Uncertainties in Radiative Transfer Computations – Consequences on the Meris Level-2 Products –

Operational MERIS level-2 processing uses auxiliary data generated by two radiative transfer tools. These two codes simulate upwelling radiances using different approaches based on the matrix-operator method (FUB), the discrete ordinate method and the successive orders technique (LISE). Intervalidation of these two radiative transfer tools was performed in order to implement them in the MERIS level-2 processing. The parameterization of the water vapor absorption defined for each of these two codes leads to a well agreement not only for the MERIS bands with residual absorption but also in the MERIS band centred at 900 nm which is used for the water vapor retrieval. As for the strong oxygen absorption at the 760.625 nm MERIS wavelength, its parameterization varies between the two codes. Nevertheless, the systematic biases in the two codes will be removed thanks to the use of a differential method between two MERIS adjacent bands. For the oxygen absorption at 760.625 nm, a more exhaustive study needs to be achieved. An extensive exercise was also conducted for cases without gaseous absorption. The scattering processes both by the molecules and the aerosols were retrieved within few tenths of a percent. Nevertheless, some substantial discrepancies occurred if the polarization is not taken into account mainly in the Rayleigh scattering computations. Over land, errors on the aerosol optical thickness reach up to 30 percent in some geometries as observed in the SeaWiFS images. Over ocean, the impact of the code inaccuracy in the retrieval of water-leaving radiances is large, up to 50% in relative difference. Applying the OC2 algorithm, the effect on the retrieval of the chlorophyll concentration is less than 10%. 1. RADIATIVE TRANSFER CODES (RTC) For the molecular scattering, the Rayleigh theory applies to the response of the electrical dipole to an incident electromagnetic field. Application of the Rayleigh theory to the atmosphere results on the choice of the pressuretemperature profile (here the mid-latitude summer (MLS)) and on the molecular anisotropy. The Mie's theory gives the inherent optical properties (i.e., extinction and scattering coefficients and phase function) of a spherical aerosol. In order to well describe the selective gaseous absorption in the radiative transfer simulations, the transmission functions have been computed from absorption line parameters and averaged over the spectral response of MERIS. The absorption line parameters of the relevant gases are extracted from the HITRAN-2000 database [1] and the atmospheric absorption is calculated from the line-by-line (LBL) model with a high spectral resolution. Considering an isotropic medium illuminated by an unpolarized radiation beam, i.e., with the Stokes parameters ) 0 , 0 , 0 , ( o E and with symmetrical boundary conditions with respect to the incident plane, the Stokes vector components can then be developed in azimuth using a Fourier series expansion [2]. By expanding the scattering phase function in the same way, the RTE splits up into a set of independent equations from which the number of equations depends on the length of the scattering phase function development (i.e., the number of terms used in the Legendre polynomial expansion). Each equation can then be solved independently in the Fourier space. The RTC from the LISE institute is based on two distinctive radiative transfer tools. The first one, namely GAME (Global Absorption ModEl), computes the radiative transfer within an absorbing atmosphere whatever the sky conditions (i.e., clearor cloudly-sky) with an approach based on the discrete ordinate method (DOM) [3]. The second one resolves the RTE for a purely scattering medium using the SO method. __________________________________________________________________________________________________________ Proc. of Envisat Validation Workshop, Frascati, Italy, 9 – 13 December 2002 (ESA SP-531, August 2003) The GAME code, which accounts for the coupling between scattering and gaseous absorption, is specially devoted to compute radiative tranfer within clouds and absorbing atmospheres. Overlapping atmospheric gaseous absorption and non-grey absorption for multiple scattering media are accurately treated with the correlated k-distribution method including continuum absorption. Transmission function ) (u T for an absorber amount u , can then be approximated by the exponential series function transmission (ESFT). However, this approximation can lead to non-negligible errors for strong absorption bands, especially in the case of the oxygen absorption for which absorption in the high atmosphere is important. Coefficients of the ESFT are calculated from a LBL model, according to the correlated k-distribution. For each term of the ESFT, the RTE is then solved which yields to the optical thickness of the corresponding atmospheric gas. Numerous tests stressed that seven exponential terms give a good compromise between accuracy and computer time. Vertical profile (33 layers) of the extinction optical thickness, single scattering albedo and terms of the Legendre polynomial decomposition derived from the optical properties of atmospheric constituents (aerosols, clouds and molecules) and gaseous absorption coefficients (H2O, O2, O3) are used as inputs to the DOM. Note that the polarization processes are not included in this code. The RTC / SO [4] solves the RTE for a purely scattering atmosphere. Over land, the atmosphere is bounded by a Lambertian ground reflector with clear-sky conditions. Over ocean, the boundary condition is a rough sea surface with wave slope orientations governed by a distribution function ) , ( n n f φ μ depending on the wind speed [5]. Radiances are then computed for each term m of the Fourier series. Scattering phase function of the aerosols are approximated by 80 terms in the Legendre expansion. Integration with respect to the cosine of the incident zenith angle ' μ is performed using a Gaussian quadrature with 48 angles and the computation with respect to the optical thickness δ is accomplished by dividing the atmosphere into 33 layers. For a given term m of the Fourier series, the iterative process over the SO of scattering ) 1 ( > n is broken as soon as the SO of scattering converges into a geometrical series. If this is the case, the iterations are then stopped and the tail of the geometrical series is added to the contribution of the scattering order for which the convergence was detected. The RTC / FUB, so-called MOMO (Matrix Operator MethOd) simulates the radiative transfer processes within a multilayered absorbing and scattering medium [6]. The matrix-operator method, also known as the doubling-adding method, relies on the assumption that the vertical structure of the medium can be approximated by an appropriate number of homogeneous layers. The development of this method is fundamentally based on the single scattering approximation. The interaction principle relates the outward directed light field at the layer boundaries of an absorbing, scattering and emitting medium linearly to the light field impinging on the layer boundaries and the radiation generated inside the layer. Distinguishing between upward and downward directed radiation, this can be expressed by reflection, transmission and source operators. The atmospheric gaseous absorption is treated with a modified k-distribution method [7]. The transmission function is approximated and included in the process of solving the RTE in multiple scattering. Thus, the interaction of multiple scattering and gaseous absorption are accurately treated. 2. RTC AND GASEOUS ABSORPTION Computations of the gaseous absorption account for the line absorption (H2O and O2) with strong spectral variations as well as the continuum absorption (H2O and O3). Outside of the strong absorption bands (i.e., in the 765.625 nm MERIS band for O2, and 900 nm for H2O), the coupling between scattering and gaseous absorption remains relatively weak. This assumption leads to express the apparent reflectance * ρ at top of the atmosphere (TOA) as the product of * na ρ , the signal ignoring the gaseous absorption, and g T the gaseous transmittivity. Correction for ozone is classical. The total water vapor content can be estimated using the radiance ratio between the 900 nm (with absorption) and 885 nm (without absorption) MERIS wavelengths. The water vapor transmittance, mostly at 708.75 nm, within a slightly contaminated spectral band can be directly evaluated using the relationship with the 900 nm to 885 nm radiance ratio. In the same way, the oxygen transmittance at the 778.75 nm MERIS wavelength is expressed versus the 760.625 nm to 753.75 nm radiance ratio. The water vapor transmittivity ) ( O H2 T was computed using the two RTCs (GAME and MOMO) in the MERIS band at 900 nm. For MERIS, this transmittivity is derived from the ratio of measurements acquired at 900 nm and 885 nm. This ratio is used to determine the water vapor content [8] [9]. A reference value for O H 2 T is given by a LBL code [7] and a comparison is reported in Fig. 1. FUB results well match the LBL computations while a systematic bias of one percent exists with GAME computations. This one percent error on O H 2 T will lead to an error of about 10 percents on the water vapor retrieval, and the MOMO code will be then used in the water vapor retrieval. Similar computations were performed at 708.75 nm for which we have the largest contamination by H2O and the results reported in Fig. 2 indicate a good retrieval of O H 2 T by the two RTCs. As previously, the oxygen transmittivity ) ( 2 O T was also computed with GAME and MOMO but in the MERIS band at 760.625 nm. For MERIS, this transmittivity is obtained from the ratio of measurements acquired at 760.625 nm and 753.75 nm. A reference value for 2 O T is given by a LBL code and a comparison is reported in Fig. 3. At 760.625 nm, the absorption by O2 is strong and the two codes do not agree well with the LBL computations. FUB results well describe the variations with the airmass but with a bias. The discrepancies are much higher with the GAME code. The same kind of computations was done at 778.75 nm for which we have a residual contamination by O2 and the results indicated a good retrieval of 2 O T with the RTC/FUB (MOMO) and an acceptable retrieval with the RTC/LISE (GAME) (Fig. 4). 2 3 4 5 6 7 0.5 0.55 0.6 0.65 0.7 0.75 Airmass H 2O T ra ns m itt ivi ty a t 9 00 nm LBL RTC/FUB RTC/GAME 2 3 4 5 6 7 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 Airmass H 2O T ra ns m itt ivi ty a t 7 08 .7 5n m LBL RTC/FUB RTC/GAME Fig. 1: vapor transmittivity at 900 nm computed with the LBL, RTC / FUB (MOMO) and RTC / LISE (GAME), versus the airmass. Fig. 2: Same legend as Fig. 1, but at the 708.75 nm MERIS wavelength. 2 3 4 5 6 7 0.2 0.25 0.3 0.35 0.4 0.45 Airmass O 2 Tr an sm itt ivi ty a t 7 60 .6 25 nm LBL RTC/FUB RTC/GAME 2 3 4 5 6 7 0.99 0.991 0.992 0.993 0.994 0.995 0.996 0.997 0.998 0.999 1 Airmass O 2 Tr an sm itt ivi ty a t 7 78 .7 5n m LBL RTC/FUB RTC/GAME Fig. 3: Oxygen transmittivity at 760.625 nm computed with the LBL, RTC / FUB (MOMO) and RTC / LISE (GAME), versus the airmass. Fig. 4: Same legend as Fig. 3, but at the 778.75 nm MERIS wavelength. The oxygen absorption band at 750 nm is commonly used to extract the cloud top pressure [10] as well as the land surface pressure [11]. Alternatively, a test on this pressure determination allows to identify clouds. Without scattering, the oxygen absorption is accurately computed with a LBL code [12] accounting for the spectral position and shape of the lines from a spectroscopic database. From satellite data, the oxygen absorption can be achieved from the radiance ratio using a spectral band in the oxygen A-band and a close non-absorbing channel. This radiance ratio is corrected from the coupling between scattering and oxygen absorption using a corrective factor C [11]. Computation of the corrective factor C to account for the coupling between scattering and absorption are depicted on Fig. 5. Discrepancies on C are around one percent within the MERIS field of view of 40 degrees when the sun is lower which has an impact of about 15 hla on the pressure retrieval. Fig. 5: Left side: 3D representation (x-scale: view zenith angle VZA, y-scale: view azimuth angle from the solar plane RAA) of the correction factor (C ) derived from the RTC / FUB (MOMO). computed at sea level for a continental aerosol model ) 33 . 1 ; 4 ( = = m α and a visibility of 23 km. Computations are performed over a reflective surface ) 2 . 0 ( = G ρ and for a SZA of 30°. Right side: relative differences on C values derived from RTC / FUB (MOMO) and RTC /LISE (GAME). 3. RTC AND SCATTERING The results stress a well agreement between the two Mie's computations, with a few 10 percent of mean relative differences which remain insignificant. The phase function computations will result in any lost of sensitivity for the atmospheric correction algorithm over land. This directly concerns the calculation of the primary scattering in the aerosol layer but also the different coupling reflectances for the dense dark vegetation (DDV). The SO code allows to compute the radiance accounting for the polarization. The available version of the MOMO code does not include the polarization processes. A comparison of the two RTCs is then conducted without accounting for the polarization in the SO code. We ran 4 test cases with a pure molecular layer (a) and with a pure aerosol layer according to the structure of the land algorithm [13] (b) plus two TOA computations with a realistic atmosphere over a dark surface (c) and over a bright surface (d). As indicated in Table 1, the two RTCs compared within better that one percent or so when only the scattering processes are involved, the code inaccuracy will have no impact on the MERIS products. Table 1: Absolute (x 10) and relative differences (max. and mean values expressed as percent) between LISE and FUB derived upwelling radiances simulated over 4 land test cases at 412.5 nm: (1) pure molecular atmosphere 314 . 0 = R δ , (2) pure aerosol atmosphere with a continental model 23 . 0 ) 550 ( = nm a δ , (3) molecules and the previous aerosol model. The first 3 test cases are over a black surface while the case (d) corresponds to case (c) with a reflective surface ) 5 . 0 ( = G ρ . The sets of three results correspond to the solar zenith angles of 7°, 30° and 70°, respectively. Case max L ∆ max L L / ∆ mean L ∆ mean L L / ∆ (a) 3, 7, 5 0.09, 0.24, 0.28 2, 2, 2 0.05, 0.07, 0.08 (b) 10, 10, 32 1.37, 1.52, 1.55 4, 4, 6 0.60, 0.66, 0.90 (c) 14, 25, 77 0.33, 0.57, 1.61 9, 10, 13 0.21, 0.25, 0.41 (d) 28, 26, 73 0.16, 0.18, 0.95 19, 16, 15 0.11, 0.11, 0.24 A set of 8 black ocean test cases was tested over two levels of surface roughness (i.e., induced by a wind speed of 3 m/s and 7.2 m/s) with different atmospheric conditions (Table 2). Both in the two RTCs (FUB and LISE), the atmosphere can be divided in three aerosol sub-layers (one aerosol type per sub-layer). For this RTC / intervalidation exercise, three aerosols models were defined: (1) the 'aerosol model 1' built with one maritime aerosol model, (2) the 'aerosol model 2' which is composed of three models, i.e., a maritime model (log-normal size distribution) for the tropospheric layer, a continental model for the second layer and a stratospheric model for the upper layer, and (3) the 'aerosol model 3' which is similar to the 'aerosol model 2' but with a maritime model (Junge size distribution) for the tropospheric layer. The visibility is set to 23 km. The simulations are performed for three solar zenith angles (SZA): 7°, 30° and 70°. The radiances were compared for a view zenith angle (VZA) ranging from 0° to 60° and a relative azimuth angle (RAA) ranging from 0° (forward direction) to 180° (backward direction). Table 2: Black ocean test cases. Case 1 2 3 4 5 6 7 8 Wavelength (nm) 412 412 412 708 412 412 708 412 Aerosol scattering No 2 1 1 1 3 3 3 Wind speed (m.s) 7.2 7.2 7.2 7.2 3.0 7.2 7.2 3.0 The maximum and mean absolute (resp., relative) differences for all the 8 black ocean test cases are given, for three SZA values, in Table 3. Table 3: Absolute and relative differences (max. and mean values) between LISE and FUB upwelling radiances over black ocean test cases for the three SZAs ) ( s θ . The test case 1 shows that the effect of the air-sea interface on Fresnel reflection and refraction are well implemented both in the two RTCs. In fact, both for the three SZAs, we observed that the bi-directional radiance distributions (BRD) of FUB and BRD of LISE present the same wavy variation with a local peak radiance in the specular direction for the two smallest SZAs [14]. The slight discrepancies (i.e., max L L | / | ∆ < 1.2%) observed between FUB and LISE upwelling radiances at TOA stress this excellent agreement between the two RTCs. As expected, above a realistic oceanic atmosphere (Rayleigh + aerosol scattering), the VIS BRDs simulated by the two RTCs (FUB and LISE) over the same roughness black sea surface are always in agreement, whatever the number of aerosol layers used as well as the optical properties selected (test cases 2, 3 and 6). Mean relative deviations between the simulated BRDs over these black ocean test cases are lower than 0.3%, 0.8% and 1.2% at 7, 30 and 70 SZA respectively. However, in the NIR region with the same atmospheric conditions over this surface roughness (test cases 4 and 7), these mean deviations are more significant, i.e., up to 0.7%, 1.5% and 3.3% at 7, 30 and 70 SZA respectively. As the contribution of the Rayleigh scattering is lower in NIR, differences are explained with the number of terms used in the Legendre expansion of the aerosol phase function. In the RTC / LISE, this number of terms (80 terms) is not always enough to model the phase function, in particular in the forward scattering peak and in the backscatter region. The oscillations observed on the phase function appear then in the simulated BRD pattern [14]. For an oceanic atmosphere over a slightly roughened-wind black sea surface (i.e., using a wind speed of 3 m/s), the BRDs simulated are also in agreement with mean relative deviations lower than 0.2%, 0.7% and 1.1% (test cases 5 and 8). This result confirms the well introduction of reflection processes at the air-sea interface in the two RTCs independently on the level of surface roughness. The impact of the polarization on the total radiance occurs after the primary scattering mainly for highly polarized scatters. This is the case at shorter wavelengths and for the Rayleigh scattering. In order to investigate the impact of the polarization, we ran the SO code with and without polarization. The aerosol remote sensing algorithm uses two bands in the blue at 412.5 nm and 442.5 nm. Computations of the upwelling reflectances at 442.5 nm in the principal plane for three solar angles give an averaged value of 0.2 for the reflectance. The Rayleigh scattering dominates in this spectral range and an accurate computation has to be achieved. Fig. 6 reports the differences on the upward reflectances at 442.5 nm computed with and without polarization. The y-axis range is included within ± 1 percent in reflectance with by far above what we can accept. If we ignore the polarization in our computations, the radiances are underestimated in the backward directions and overestimated at 90 degrees scattering angles. The effect of the polarization in the computation of the aerosol path radiance is illustrated in Fig. 7. Only the aerosol layer was considered with a continental model (Junge size distribution) and a meteorological visibility of 23 km. Because the aerosols less polarize than the molecules, the errors on the total reflectance when neglecting the polarization become negligible. Fig. 6: Difference in reflectances at the 442.5 nm MERIS wavelength computed at the top of a pure molecular atmosphere with the RTC / SO including or not the polarization processes Fig. 7: Same as Fig. 6 but at the top of a continental aerosol layer (visibility of 23 km). Simulations are completed in the principal plane for three SZAs: 28° (+), 58° (*) and 73° (o). Backscattering directions correspond to negative VZA values 4. IMPACT ON THE MERIS LEVEL-2 LAND PRODUCTS Assuming the same optical thickness and the same asymmetry factor, we saw that, see test case (a) from Table 1, the Rayleigh scattering is well computed by the two RTCs (FUB and LISE). The major impact results on the introduction of the polarization or not to compute the LUTs and mostly occurs in the blue. To illustrate this effect, we applied the MERIS-like algorithm to a SeaWiFS image acquired over the South West of France on September 03, 1999 at 12h40 UT. Fig. 8 displays the histogram of the SeaWiFS reflectances at 412.5 nm after a Rayleigh correction accounting or not for the polarization. The solar angle was 40 degrees for a scattering angle of 154 degrees for the main DDV area. For this specific geometry, the Rayleigh reflectance is higher by about one percent when we account for the polarization. The first stage is to select DDV pixels using the atmospherically resistant vegetation index (ARVI). The latter is obtained after a Rayleigh correction. Because the Rayleigh correction is different between accounting or not for the polarization, the DDV selection differs. Additional DDV pixels appear when ignoring the polarization which also illustrates the sensitivity of the ARVI threshold on the selection. The AOT a δ at 443 nm and 650 nm are then derived from the DDV pixels. At 443 nm, the a δ determination directly results from the difference in the Rayleigh correction and reaches here 0.05 (Fig. 9). Clearly, the aerosol products are significantly biased when the polarization is not accounted for. In the red region at 670 nm, the Rayleigh reflectance computation is no longer subject to the polarization. The spectral dependency of a δ is then stronger when the polarization is accounted for. That impacts on the aerosol model selection. We also investigated the error on the surface reflectances. The atmospheric path radiance is correctly removed in the two cases, because this correction directly results on the DDV reflectance retrieval. The differences between with and without polarization originate from the computation of the aerosol transmittances. Fig. 8: Histogram of the aerosol reflectance (Rayleigh effects removed) at 412.5 nm and 670 nm from the SeaWiFS scene, with (full line) and without (dashed line) polarization Fig. 9: Histogram of the aerosol optical thickness at 443 nm estimated over SeaWiFS scene, with (full line) and witout (small dashed line) polarization 5. IMPACT OF THE POLARIZATION ON THE OCEAN LEVEL-2 PRODUCTS The best way to look at the sensitivity of the ocean colour product to the effect of the polarization of the light field is to generate a second set of LUTs accounting for the polarization in the RTC simulations. It was not possible for this paper because of consuming computer time. For a preliminary study, we developed an approach based on the 6δ like formulation of the signal. The TOA radiance toa L3 computed with the polarization for the atmosphere-ocean system for which the water leaving radiance retrieved w L is known. Then toa L is inverted with atmospheric functions computed without polarization to get an other value unpol w L of the water leaving radiance. The difference between retrieved w L and unpol w L represents the impact of the consideration of the polarization on the level-2 product. We quantify the impact of the chlorophyll ) (Chl concentration product using the OC2 algorithms [15]. A set of 4 coupled 'Atmosphere-Ocean' cases was tested (Table 4). We define two aerosol models, a continental one and a maritime one. The simulations are performed for two SZAs (30°, 60°) representing respectively summer and winter conditions. The radiances were simulated for three VZAs (0°, 18.5°, 41.4°) and for 4 RAAs (0°, 45°, 135°, 180°). Table 4: Test cases for a coupled 'Atmosphere-Ocean' system. Absolute and relative differences (max. and mean values) between retrieved w L and unpol w L The differences between retrieved w L and unpol w L are large whatever the test cases, up to 47%. In the NIR region, this difference decreases with the Rayleigh scattering contribution (test cases 6, 9 and 12). The polarization must be taken into account in the LUT generation for the atmospheric correction. The relative maximum difference between the 2 Chl concentrations retrieved from retrieved w L and from unpol w L using OC2 algorithm is 10% (Table 5). In fact, the error on the Chl concentration is 0.2 μg/l for a Chl concentration of 2 μg/l and 0.05 μg/l for a Chl concentration of 0.5 μg/l. The OC2 algorithm based on the ratio of the reflectance at 490 nm and 560 nm seems to be robust to the inaccuracy in the simulation of the atmospheric radiance and path without polarization. Table 5: Absolute and relative differences (max. & mean values) between Chl concentration retrieval from retrieved w L and Chl concentration retrieval from unpol w L using OC2 algorithm.