Speeding up the aerosol optical thickness retrieval using analytical solutions of radiative transfer theory

We present here the aerosol retrieval technique FAR that uses radiative transfer computations in the process of retrieval rather than look-up tables (LUT). This approach provides operational satellite data processing due to the use of the accurate and extremely fast radiative transfer code RAY previously developed by authors along with approximate analytical solutions of the radiative transfer theory. The model of the stratified atmosphere is taken as two coupled layers. Both layers include aerosol scattering and absorption, molecular scattering and gas absorption. The atmosphere parameters are assumed to change from pixel to pixel in the lower atmosphere layer, but the upper stratified layer of the atmosphere over 2–3 km is supposed to be horizontally homogenous for the frame under retrieval. The model of the land spectral albedo is taken as a weighted sum of two a priory chosen basic spectra. The aerosol optical thickness (AOT), Angstr öm exponent and the weight in the land spectral albedo are optimized in the iteration process using the least-squares technique with fast computations of the derivatives of radiative characteristics with respect to retrieved values. The aerosol model and, hence, the aerosol phase function and single scattering albedo, is predefined and does not change in the iteration process. The presented version of FAR is adjusted to process the MERIS data. But it is important that the developed technique can be adapted for processing data of various satellite instruments (including any spectral multi-angle polarizationsensitive sensors). The use of approximate analytical radiative transfer solutions considerably speeds up data processing but may lead to about 15–20% increase of AOT retrieval errors. This approach is advantageous when just the satellite data processing time rather than high accuracy of the AOT retrieval is cruCorrespondence to: I. L. Katsev (katsev@light.basnet.by) cial. A good example is monitoring the trans-boundary transfer of aerosol impurities, particularly in the case of emergencies such as volcano eruptions, or various industrial disasters. Beside, two important problems that determine the accuracy of the AOT retrieval are considered. The first one is the effect of the preliminary choice of the aerosol model, particularly for the retrieval from satellite instruments providing only spectral data (MERIS, MODIS). The second problem is the influence of clouds in adjacent pixels. As for our knowledge, this problem has not been given the required attention up to now and it should be properly accounted for in the AOT retrieval algorithms.


Introduction
For many years, all main techniques to retrieve a spectral aerosol optical thickness (AOT) from satellite data include the radiative transfer (RT) in the data processing using lookup tables (LUT).The main and evident advantage of this technique is time saving.But increasing information provided by satellite sensors (multi-angle, multi-spectral data, polarization measurements) opens the possibility of using various statistical optimizations in satellite data processing.With this in mind, the use of RT calculations in the satellite data processing has great potential.This is why, a few years ago, we started to develop the aerosol retrieval techniques that use radiative transfer computations in the process of retrieval rather than LUT.The first operating algorithm of this type, the ART code for processing MERIS spectral data, was presented by Katsev et al. (2009).Other codes, that use RT calculations for satellite data processing, have appeared lately as well (Kokhanovsky et al., 2010a).
But such retrieval techniques can be applied operationally only if the accurate and extremely fast radiative transfer procedures are used.Previously, we have developed the RAY code (Chaikovskaya et al., 1999;Tynes et al., I. L. Katsev et al.: Speeding up the AOT retrieval procedure using RTT 2001) for simulations of radiative transfer in the atmosphereunderlying surface system with regard to polarization.This code meets accuracy and speed requirements.The short description of RAY was given in Katsev et al. (2009).Recently its accuracy, along with high calculation speed, was confirmed in Kokhanovsky et al. (2010b).For instance, the relative error of RAY for the first Stokes vector component is estimated to be about 0.003% for Rayleigh scattering and about 0.2% for aerosol scattering at AOT = 0.3.The RAY's high processing speed allows the use of the iterative radiation transfer computations in satellite data processing for the AOT retrieval.The RAY code is a core of codes for the AOT retrieval developed in the B.I. Stepanov Institute of Physics of National Academy of Sciences of Belarus (see, for instance in Katsev et al., 2009).
The first successful implementation of this idea (the use of RT computations in the process of the AOT retrieval with the least-squares optimization) was the ART code for MERIS data processing outlined in Katsev et al. (2009).At the moment, the fundamental goal for the satellite aerosol monitoring is to achieve the accuracy of an aerosol retrieval technique that meets the requirements of climate problems.Nevertheless, there exists some very important problems where just the satellite data processing time rather than high accuracy of the AOT retrieval is crucial and the accuracy of the AOT retrieval about 10-20% is acceptable.A good example is monitoring the trans-boundary transfer of impurities, particularly in the case of emergencies such as volcano eruptions and various industrial disasters.Below, we present a new FAR (Fast Aerosol Retrieval) code that has a lot in common with the previous ART code.The main difference is in the used RT procedures.Like the ART code that can be considered as a prototype, the FAR does not use the LUT technique.Unlike the ART, it includes the approximate analytical solutions for the lowest atmosphere layer characterised by the highest spatial variability, along with the RAY computations for the atmosphere above this layer.This new approach keeps the main advantages of the ART technique, allowing a lot of convenient features, for instance, simple revising of atmosphere models and of aerosol composition, use of the least-squares method to find the AOT and Angström exponent from satellite data and so on.You will see in Sect.3, the errors of the AOT retrieval with the FAR only slightly (about 15-20%) exceed errors provided by the ART retrieval and stay less than inaccurate due to other factors including a priori choice of optical models of aerosol and surface.But retrieval with the FAR is about 100 times faster than with the ART technique.
Both codes can be used for the AOT retrieval for atmosphere over land and over water, but here for brevity we will consider only retrieval of the AOT over land.
This paper is arranged as follows: The FAR algorithm is presented in Sect. 2. The very short description of the used atmosphere and surface models of the procedures of the FAR code, same as in the ART (preparation of the in-put data, optimization with the least-squares technique), are given in Sect.2.1.Because all details can be found in Katsev et al. (2009), this section only includes information that is necessary for the understanding of the FAR technique.Section 2.2 introduces the approximate analytical solutions used in the FAR for the calculation of the radiative characteristics of the lower aerosol layer and illustrates their accuracy.The radiative interactions between layers are considered in Sect.2.3.Section 3 is devoted to the validation of the FAR code with the ART for benchmarking and by comparison with AERONET data.Besides, two important problems that determine the accuracy of the AOT retrieval are considered.The first one is the effect of the preliminary choice of the aerosol model, particularly for the retrieval from satellite instrument providing only spectral data (MERIS, MODIS) (Sect.4.1).The second problem is the influence of clouds in the adjacent pixels (Sect.5).As for our knowledge, this problem has not been given the required attention up to now and it should be properly accounted for, practically in all AOT retrieval algorithms.

Aerosol and surface models and brief description of the AOT retrieval
The developed FAR algorithm is designed to retrieve the AOT, Angström exponent and underlying spectral surface albedo from the top of atmosphere (TOA) spectral reflectance defined as where E 0 (λ) is the extraterrestrial irradiance incident normally on a given unit area at the top-of-atmosphere, I ↑ λ,µ,µ 0 ,ϕ is the measured TOA radiance, µ 0 and µ are cosines of the incidence and observation zenith angles, ϕ is a difference between the azimuth angles of the incidence and observation directions.
As it was mentioned above, the ART algorithm described in detail in Katsev et al. (2009) is a prototype of the presented FAR algorithm.One can consider this FAR algorithm as a faster version of the ART.In this section, we will give a flow chart of the FAR algorithm briefly reminding the features of the data processing that are common for both codes (ART and FAR) and referring a reader to the abovementioned publication for the detail.
In both algorithms, the model of the stratified atmosphere is taken as two coupled layers.Both layers include aerosol scattering and absorption, molecular scattering and gas absorption.The layer "1" (lower layer) of the atmosphere is a comparatively thin layer of the lower troposphere up to the height H .As a rule, the value of H is about 2-3 km.As the aerosol in this layer is characterised by maximal spatial and temporal variations, the AOT of this layer is supposed to vary from pixel to pixel and is retrieved for each pixel independently.Stratification in this layer is neglected.The performed computations showed that neglecting the stratification inside layer "1" leads to the relative error of calculations of the reflectance at the top-of-atmosphere less than 0.2% for any wavelength in the visible.
The optical characteristics (optical thickness τ , single scattering albedo ω, and phase function p(θ)) of this uniform layer "1" are calculated including contributions of molecular scattering, gas absorption and absorption and scattering by aerosol, namely where τ R , τ aer and τ g are the optical thicknesses corresponding to molecular scattering, light extinction by aerosol particles and gaseous absorption contributions, respectively.Equation (3) also includes the single scattering albedo of aerosol particles (ω aer ) and that of Rayleigh scattering processes (ω R = 1).The coefficients x n of the expansion of the phase function p(θ) into series of Legendre polynomials are Here x n,R and x n,aer are the coefficients of the expansion of the Rayleigh and aerosol phase functions, respectively, into the series of Legendre polynomials.These coefficients x n are used in the adding method procedure to calculate the radiative interaction between layers "1" and "2".Because the Rayleigh scattering phase function comprises only of the zeroth and second Legendre polynomials, we have: x n = 0.5δ n2 τ R + ω aer τ aer x n,aer τ R + ω aer τ aer (5) at n ≥ 1 and also x 0 = 1 by definition.Here δ n2 is the Kronecker symbol equal to one at n = 2 and zero, otherwise.The layer "2" (upper layer, i.e. the atmosphere above the altitude H ) includes the stratosphere and the upper and middle troposphere.Naturally, the layer "2" is characterised by the vertical stratification of the aerosol and gases concentrations, pressure and temperature profiles.In the computation procedure, this layer is considered to be composed of N homogenous sublayers with optical characteristics averaged over each sublayer.The optical characteristics of each sublayer are computed including contributions of molecular scattering, gas absorption and aerosol absorption and scattering similar to the procedure described above with Eqs.(2-5).To reduce the computation volume the radiation characteristics of the stratified layer "2" are computed one time for all pixels of the processing frame.For layer "2", the reflectances at the illumination from the top R 2 µ,µ 0 ,φ and from the bottom R * 2 µ,µ 0 ,φ as well as the transmittance at the illumination from the bottom T * 2 µ,µ 0 ,φ are computed with the fast and accurate RAY code.
In the iteration procedure for the retrieval of the AOT and Angström parameter, the model of the land spectral albedo r s (λ) is taken as a linear combination of some basic spectra of the vegetation r veg (λ) and soil r soil (λ).These basic spectra are given by von Hoyningen-Huene et al. (2003).Thus, the spectral albedo of the surface is characterised by the only parameter c.Just this surface parameter is determined in the retrieval process.The normalized differential indices serve as a zeroth approximation for the parameter c for each pixel.This normalized differential vegetation index for the "land" pixels is defined as In this paper, we will present examples of the retrieval for the middle latitudes and the land reflectance model ( 6) is used.
For other cases, the different basic functions in Eq. ( 6) might be used in FAR.
The spectral radiance at the TOA measured by any satellite optical instrument in N spectral channels could be the input to the FAR (and ART) algorithms.Because both algorithms do not use the LUT technique, the choice of the spectral channels is very flexible.Their number and particular wavelengths depend on the spectral characteristics of the used satellite optical instrument.Satellite data in 9 MERIS spectral channels specified by the wavelengths λ = 412.5,442.5,490,510,560,620,665,865and 885 nm are used.
The flow chart of the FAR code is presented in Fig. 1.The preparation of the input data includes operations described in steps 1-3 below.Steps 4-10 describe the sequence of the retrieval operations.
2. Separating the "land" (R TOA (885 nm) ≥ 0.1) and "water" (R TOA (0.885 nm) < 0.1) pixels (von Hoyningen-Huene et al., 2003) and discarding "bad" pixels in the group of "land" pixels on the basis of the NDVI values: pixels classified as "land" with NDVI < 0.1 are considered as pixels containing sub-pixels "water" and are discarded (Remer et al., 2005).scenario and the problem under consideration.In the particular case of MERIS data, we used the box with 5×5 pixels.For this box the average value of RTOA (λ) is calculated.In doing so 20% of the pixels with minimal values of R TOA (665 nm) and 30% of the pixels with maximal values of R TOA (665 nm) are discarded from the data to be averaged.The determined average value of RTOA (λ) is ascribed to the central pixel in the box.
4. Computations of the reflectance at the illumination from the top R 2 µ,µ 0 ,ϕ and from the bottom R * 2 µ,µ 0 ,ϕ as well as the transmittance at the illumination from the bottom T * 2 µ,µ 0 ,ϕ for the layer "2".These radiative transfer characteristics are computed taking into account the stratification of atmosphere and light polarization effects as a solution of the vector radiative transfer equation with the RAY code. 5. Specifying the starting values of the sought for parameters for the iteration process.The zeroth approximation can be chosen as: τ 412,0 = 0.3, α 0 = 1.3, and c 0 = NDVI.Practice with the AOT retrieval shows that the iteration process convergence is hardly sensitive to the choice of the zeroth approximation for the parameters α and τ 412 (Katsev et al., 2009).
7. Calculating the corrections τ 412, i , α i , and c i with least square technique.The corrections at the i-th step of iterations are determined from the following set of equations: where RTOA λ j are the averaged values of the TOA spectral reflectance obtained after the satellite data processing described in steps 1-3, the values R TOA,i λ j are the radiances computed at τ 412 = τ 412,i , α = α i , and c = c i .
8. Checking satisfiability of the condition If this condition is not fulfilled, the new values are calculated as and the iteration procedure (steps 6-8), depicted in Fig. 1 Note that in this version of the FAR code the aerosol model and, hence, the aerosol phase function and single scattering albedo, is predefined and does not change in the iteration process.9. Calculating the spectral optical thickness for the total atmosphere τ 2 (λ) being the optical thickness of aerosol in the stratified upper layer "2", and the Angström exponent for the total atmosphere.
10. Retrieving the corrected spectral albedo of the underlying surface r s (λ) with the retrieved dependenceτ (λ).
When the underlying surface is land, the surface albedo r s is defined from the well-known equation (Chandrasekhar, 1960): where R a µ,µ 0 ,ϕ is the TOA reflectance of the whole atmosphere above a black underlying surface, t a µ 0 is the transmission coefficient of the whole atmosphere, r s , the surface albedo, r * sa , the spherical albedo of the atmosphere at illumination from an atmosphere bottom upwards.
To conclude this section, let us make two following comments.
Firstly, Eq. ( 6) for the surface spectral albedo is used only in the iteration process.Finally the spectral albedo is retrieved at step 10 with the spectral AOT obtained with the iteration procedure.
Secondly, regarding the bidirectional reflectance of the underlying surface, in fact, in the iteration procedure to retrieve the AOT and Angström exponent the following equation is used for a particular pixel with fixed parameters of the geometry µ,µ 0 ,φ : Equation ( 15) includes some effective reflection coefficient r ef s instead of the albedo r s .This value depends on BRDF, the incidence and observation angles and on the ratio of diffuse to direct components of the solar light after propagating through the atmosphere.This means that the value r ef s also depends on the atmosphere parameters, which determine this ratio.Just an effective surface reflection coefficient r ef s for the particular observation geometry that corresponds to the particular pixel is included in the iteration process.It means that BRDF of the underlying surface (non Lambertian surface reflectance) does not lead to the additional AOT errors in comparison with those in the case of a Lambertian underlying surface.Note that the spherical albedo r s is supposed to stay in the denominator of Eq. ( 15) (as in Eq. ( 14) because at the multiple re-reflections between the surface and atmosphere the surface is illuminated by the diffuse light.But the value r s r * sa 1 and an error due to change the value of r s on r ef s has a little effect on the accuracy of R TOA µ,µ 0 ,φ calculation.Regarding the above note, it can be asserted that FAR algorithm allows for the surface BRDF at AOT retrieval.

Analytical radiative transfer solution used in the FAR code
To speed up satellite data processing, the FAR algorithm uses approximate analytical solutions of the radiative transfer theory to calculate the reflectance R 1 µ,µ 0 ,φ and transmission coefficient t 1 (µ) for the troposphere layer "1" considered being uniform.We emphasize that radiative transfer in the layer "2" is computed accurately with the RAY code.The main difficulty in the development of the simple and sufficiently accurate approximation is the elongation of the troposphere aerosol phase function.Actually, the solution of this problem is well-known: it is the truncation of the phase function, the simplest version being known as δ-Eddington approximation.The truncation of the phase function in small angle region does not practically affect the reflectance of the layer in the backward hemisphere.Figure 2 that shows the reflectance for the aerosol layers with optical thicknesses τ = 0.2 and 0.5 versus the observation angle illustrates this statement.This figure presents data for the layer with the Continental aerosol under the normal incidence.The calculations are performed at different truncating angles with code RAY.It is seen that even the truncation at the angle of 60 • only negligibly changes the reflectance at observation angles ϑ < 80 • .For the truncated phase function the reflectance from a layer can be calculated using any simple approximation developed for a not very elongated phase function with optical parameters τ , ω, p(θ) (the sign "∼" shows that the value is defined) instead of τ , ω and p(θ).These optical parameters are defined as: Here ω, τ and p(θ) are the single scattering albedo, optical thickness, and phase function, θ, the scattering angle, η, the truncated part of the phase function, θ * , truncating angle.Here the phase function is normalized as (1/2) π 0 p(θ)sinθdθ=1.
For the reflectance from a layer with the truncated phase function, we use simple formulas given in (Sobolev, 1975) for not very elongated phase functions.With this solution, the single scattering is computed accurately.The multiple scattering for the truncated phase function is determined as more accurate than the less elongated phase function.It is why the merging of phase function truncation with this approach allows reasonably accurate computations for layer "1".With this approximation (for the brevity let us refer to it as MSA (Modified Sobolev Approximation) in what follows) the function R 1 µ,µ 0 ,φ is calculated as Here R (1) x1 = 3 g is the first coefficient of the expansion of the truncated scattering phase functions p(θ) into a series of Legendre polynomials, g, the average cosine of the scattering angle for the truncated phase function.
As it is seen from Fig. 2, such truncation does not cause the noticeable errors in the computed reflection functions from layer "1".Simultaneously it leads to the values of x * 1 ∼ 1 for the truncated phase functions that provide the reasonable accuracy of the used approximation.Figure 3 that compares the RAY and MSA calculations of the reflectance for the layers with optical thicknesses in the range 0-1 and for various geometries of observations illustrates this statement.The analysis shows the truncation angle θ * ≈ 45 • is about optimal for the estimation of the layer reflectance.
Figure 4, that presents the estimations of the MSA errors obtained by comparison with the RAY computations for different truncation angles, illustrates this recommendation.The overview of the accuracy of the MSA solutions at different geometries of incidences and observations one can get from Figs. 5 and 6.The accurate reflectances for the error estimations were computed with the RAY code.Here the azimuth angle 180 deg at µ = µ 0 corresponds the backscattering direction.
From these data, it follows that the errors of the MSA approximation practically do not exceed 10% for τ < 0.5 and observation angles ϑ < 60 • .One can expect that errors of the optical thickness retrieval due to the use of this approximation are less than 10% as well.The largest value of errors occur near the directions of the mirror reflection at low Sun position (Sun and observation angles more than 60 • , azimuth angle 0 • ).But this geometry practically is not used in satellite observations.Note that the black areas in Figs. 5 and 6, where we do not even estimate the errors of the MSA, include only these regions.

Radiative interaction between the atmosphere layers in the FAR code
The interaction between layers is included accurately using the layer adding method (Lenoble, 1985).The reflection and transmission functions of layer "2" are computed accurately with the RAY code taking into account the atmosphere stratification and light polarization effects, the reflectance from layer "1" is determined within the MSA by Eqs.(16-20) without accounting for polarization.As it was shown in Katsev et al. ( 2009), the errors due to neglecting polarization effects for layer "1" in the most typical situations do not exceed 1%.
The transmission coefficient t a µ 0 of the atmosphere (see Eq. 14) defined with regard to multiple re-reflections between layers "1" and "2" can be given by the following approximate relation: where and the transmission coefficient of layer "1" can be taken as where Figure 7 illustrates the accuracy of this approximation presenting the comparison between the angular dependencies of the transmission coefficients t 1 µ 0 computed with Eqs.(25-26) and with the RAY code.Data for the different values of the optical thickness τ 1 of layer "1" and for the aerosol models Continental and Water-soluble are presented.It is seen that at τ 1 < 0.5 and ϑ < 70 • the error of the approximation (25) is less than 5%.Now let us consider the accuracy of Eqs.(24-25) that define the transmission coefficient t a µ 0 of the atmosphere.As it is seen from Fig. 8, these simple solutions provide the reasonable accuracy for the computations of the atmosphere transmission coefficients.Thus, with the obtained approximations we got the reasonable accuracy for the atmosphere TOA reflectance R a µ,µ 0 ,ϕ in Eq. ( 14) that gives the TOA reflectance R TOA µ,µ 0 ,ϕ from the atmosphere-underlying surface system.Equation ( 14) includes one more value that requires the attention.It is the spherical albedo of the atmosphere at illumination of atmosphere from bottom upwards r * sa .As it is seen from ( 14) this value gives only a small correction to the value of the radiance R TOA µ,µ 0 ,φ and it is enough to get an approximate estimation of this value.Because the layer "1" is considered as uniform, we can approximately assume that Here r s,1 is the reflection coefficient for layer "1" under the diffuse illumination, r * s,2 is the reflection coefficient for the layer "2" under the diffuse illumination from below.For diffuse illumination let us use the following approximations (Zege et al., 1991): for layer "1" and for the upper layer "2" with the governing Rayleigh scattering.Formulas ( 28) and ( 29) are accurate for the case of quasi-diffusion illumination of the homogeneous layer (Zege et al., 1991).Nevertheless, we use them for diffusion illumination because, we deal with comparatively weekly absorbing layers (when the quasi-diffusion radiance distribution is close to the diffusion one).Besides, as it was mentioned above, we need to get only an approximate estimation of the value r * sa by Eqs.(27)(28)(29).

Comparisons of the FAR retrievals with the ART and AERONET data
The ART retrieval was previously carefully checked by comparing with the retrieval results of other known algorithms and with AERONET data (Katsev et al., 2009).In so doing, we used the data of the inter-comparison described in Kokhanovsky et al. (2007).In this work, the performance of different algorithms was tested for a site in Europe, where multiple and near-simultaneous satellite data were available.As many as ten different algorithms for the AOT retrieval that used data taken by six satellite optical instruments currently operated in space were compared.The explored site included the cloudless ground scene in central Europe (mainly, Germany) on 13 October 2006 (10:00 UTC).The latitude range was 49 • N-53 • N and the longitude range was 7   and with other algorithms can be found in (Katsev et al., 2009).They showed that the ART data are in a good agreement with results of MISR, MODIS, MERIS BAER, and MERIS ESA retrieval algorithms.
Several AERONET instruments operated in the area at the time of satellite measurements.We also compared the ART retrieval with AERONET data.Figure 9, where the ART retrieval data are plotted versus AERONET data for λ = 550 nm, demonstrates satisfactory agreement.AERONET AOT data at 550 nm are obtained by interpolation between λ = 440 nm and λ = 670 nm.The ART values of the AOT at λ = 412.5 nm and the Angström exponent are directly retrieved, the values of AOTs at λ = 550 nm are calculated through these values.
Thus, the accuracy of the ART, where the code RAY is used to compute the radiative transfer through atmosphere, is more than satisfactory at least for the available body of satellite information (the MERIS retrieval was considered).Now we can use the ART retrieval for the FAR benchmarking.Actually, in so doing we primarily check how the use of the approximate analytical description of the radiative transfer in layer "1" affects the accuracy of the AOT retrieval.Figure 10 presents satellite RGB images obtained from European ENVISAT platform on 11 October 2005 (left) and 1 May 2006 (right).The territories of Belarus, Baltic countries and Poland are seen in the left picture.It is the anticyclone conditions, the weather is fine, atmosphere is clear.The right image shows the territories of Belarus and west Russia.In this particular day (1 May 2006) the fire smocks caused high atmosphere turbidity.Figure 11 presents the AOT distribution over presented in Fig. 10 regions retrieved with the FAR code.
Figures 12 and 13 show very good correlation between values of the AOT retrieved with the algorithms ART and FAR using MERIS on ENVISAT data for the area specified in Fig. 10.But for small values of the AOT (τ < 0.5) the FAR retrieval provides somewhat overestimated AOT values.On the contrary, at τ > 0.5 the FAR retrieval leads to small underestimations of the retrieved AOT values.It is seen evidently from Figs. 14 and 15 where the histograms of the distributions of the optical thickness at the wavelengths 412.5 nm and 550 nm are shown.The reason of these deviations is evident: it is the approximate nature of analytical formulas used in the FAR for the description of the radiative transfer in layer "1".The accuracy of this approximation depends on the optical thickness of layer "1" and the observation geometry.Nevertheless, values of the AOT retrieved with the FAR technique differ from the correspondent the ART retrieved AOT values no more than 15-20%.This deviation is the payment for the retrieval fastness.
Let us now consider the correlation of both the FAR and ART retrieval with AERONET measurements.Figure 16 shows this comparison using data of the AERONET stations in Belsk (Poland) and Zvenigorod (near Moscow) during periods from March to September in 2008 and 2009.The ART and FAR results are pretty close and agree with AERONET data.At small AOT values of the FAR retrieval show the same bias as was observed earlier.Figure 16 depicts some bias even for the ART retrieval.
To conclude this section, let us underline the FAR potentialities as a tool for the atmosphere correction and retrieval of the land spectral albedo.The correlations between spectral albedo retrieved with the ART and FAR codes for the wavelengths 412.5 nm and 560 nm are depicted in Figs. 17 and 18. Practically both codes retrieve the same values of spectral albedo.This conclusion seems to be important because it emphasizes the potentialities of the FAR code for atmosphere correction of satellite images.It is important that the FAR noticeably speeds up data processing in comparison with the ART.Processing time for 10 6 pixels at ordinary PC with the dual-core processor (Intel Core 2 Duo E6600, 2.4 GHz) is about 3 h with the ART and 2 min with the FAR.Thus, the FAR is about 100 times faster than the ART.

Importance of the choice of the aerosol model
As mentioned above, a very important issue is a priori choice of the aerosol model of layer "1" particularly for the satellite optical sensors providing only spectral information.This model assumes not only the value of the single scattering albedo, but also what is the most important, the phase function as well.We have already touched this problem in Katsev et al. (2009) but it needs more detail consideration and clear understanding.In this section, we will consider how a priori choice of aerosol model effects the retrieved AOT and spectral surface signatures.

Influence of the aerosol model on the retrieved AOT
Preliminary choice of aerosol model is a necessary step in any AOT retrieval algorithm.The correction of this model in the process of the retrieval depends on the information, provided by satellite instrument, a priori information about local aerosol and the used retrieval technique.In all cases the preliminary choice of aerosol model effects the retrieval.This influence is more important as the less information is provided by the satellite sensor.Hence, it is particularly important for the satellite sensors providing only spectral information (MERIS, MODIS).In this section, we will consider just this case.
Let us estimate the influence of the predefined phase function on the retrieved AOT values.Figure 19 demonstrates the phase functions of the aerosol models Continental, Water-Soluble and Oceanic at λ = 550 nm and the phase function that was measured earlier in Germany (von Hoyningen-Huene et al., 2003) (we will refer to it as the experimental phase function).This phase function is used, for instance, in the retrieval code BAER (von Hoyningen-Huene et al., 2003).As seen, in the angle range ∼ 100 ÷ 150deg mainly used for the AOT retrieval the values of this experimental phase function are about two times larger than the values of the Continental phase function.
Values of the AOT retrieved with the ART algorithm with these two different aerosol models (Continental and experimental) are compared in Fig. 20.In what follows we will mainly use the retrieval with the ART technique as more accurate.The experimental model is rather simple: the single scattering albedo is equal to 0.9 (as for the Continental model) and the phase function is considered independent from the wavelength.The correlation between these data is rather high, but the AOT values retrieved with the experimental phase function are 1.5 times smaller.
This difference is easy to understand.Let us consider the simplest case of a small AOT, when the single scattering approximation can be used.In this case the reflectance of the aerosol layer "1" is proportional to the product of the single scattering albedo ω, phase function p(θ) and the AOT τ , i.e.R 2 µ,µ 0 ,ϕ ∼ ωp(θ)τ. (31) It means that, for example, a twofold increase of the value of the phase function p(θ) in a considered direction leads to a twofold decrease of the retrieved values of the AOT.Thus, the choice of the model of the troposphere aerosol in layer "1" can considerably affect the retrieved AOT values.
In Kokhanovsky et al. (2010a) the following computer experiment is described.The spectral response of the atmosphere with the black underlying surface was computed with the SCIATRAN code for a chosen atmosphere model and a few values of the AOT.The computed spectral reflectance simulated signals registered by satellite sensors.The simplest case of the black underling surface was considered.Then, the spectral AOT was retrieved with different retrieval techniques using these simulated spectral radiance coefficients.The atmosphere model and particularly aerosol microphysical (and optical) properties used in simulations of the radiances at the atmosphere top were kept unknown for   the retrievers.The results that the AOT retrieved by different groups with different techniques are given in Kokhanovsky et al. (2010a).
The results of such retrieval for three aerosols models (Continental, Water-soluble, Oceanic) are depicted in Fig. 21, dependent on the reference AOT.The phase functions for these aerosol models at λ = 550 nm are given in Fig. 19.As seen, these phase functions are very different, particularly Oceanic and Water-soluble in the angular range 100 • -150 • typical for the satellite remote sensing.It is clear that just this difference provides such a discrepancy of the retrieved AOT values.The retrieval with aerosol Oceanic is very close to the reference values in the computer experiment.Now let us consider the results of the following computer experiment.The spectral reflectance at the top-ofatmosphere were calculated with the accurate code RAY for the atmosphere models used for the retrieval for two cases (Oceanic and Water-soluble aerosols) with correspondent retrieved values of the AOT.The results are presented in Fig. 22.Let us note that the atmosphere model (that includes the models of the molecular atmosphere, gases, stratification of all components) used for the AOT retrieval could essentially differ from the model used for the simulation of the satellite data.But as seen from Fig. 22, the calculated spectral reflectance for both used aerosol models do not differ and coincide with the reflectance at the top-of-atmosphere, simulated with an unknown atmosphere model.
From the results of the described computer experiments it follows: 1.The spectral dependence of the reflectance measured at the TOA in the visible and near IR does not contain enough information for the unique choice of the aerosol model even when spectral reflectance of the underlying surface is given.The spectral measurements in the range 400-900 nm can not be used to recognize the aerosol type (at least for the numerical experiment described by Kokhanovsky et al., 2010a).In this paper, it was mentioned that there is a chance that the aerosol type can be better constrained if shortwave IR measurements are used.
2. In the best case the spectral dependence of the TOA reflectance can allow one to determine some parameters of the aerosol model if this model is chosen a priori.
3. A priori choice of the aerosol model affects noticeably the retrieved value of the AOT.

Influence of the aerosol model on the retrieved surface albedo
One of the most important problems of satellite remotesensing is the retrieval of the true spectral albedo r s (λ) of the underlying surfaces (atmospheric correction).After retrieving τ (λ) the spectral albedo r s (λ) of the underlying Lambertian surface can be simply calculated from Eq. ( 14) if the surface is land: Let us study how the choice of the aerosol model, particularly the aerosol phase function, influences the spectral surface albedo retrieval.The comparison of the spectral surface albedo for two wavelengths retrieved with Continental and experimental phase functions is shown in Fig. 23.The correlation of these data is very high.Hence, the following important conclusions from this study: the retrieved spectra of the surface albedo are comparatively stable whereas the retrieved values of the AOT are sensitive to the choice of the aerosol model.This insensitivity of the retrieved spectra of the surface albedo to variations of the aerosol phase function p(θ) follows immediately from the fact that the reflectance of layer "1" is proportional to the product of the phase function p(θ) and AOT (see Eq. 31).This statement is more accurate as the optical thickness of layer "1" is smaller.Indeed, for a thin layer a twofold increase of the value of the phase function p(θ) in the considered direction leads to a twofold decrease of the retrieved values of the AOT.The product p(θ)τ and the value of R a µ,µ 0 ,ϕ in Eq. ( 31) change only slightly.At small AOT such changes of τ lead to modest variations of the transmission coefficient t a (µ).

Estimation of influence of adjacent cloudy pixels on the retrieved AOT value
In our ART and FAR algorithms, the pixels with R TOA (560) ≥ 0.   inequality immediately discards pixels containing comparatively thick clouds.The choice of the often used threshold criterion as R TOA (560) ≥ 0.2 might lead to discarding not only cloud pixels, but pixels with highly reflective surfaces (for instance, sands, deserts).However, the choice of the threshold criterion as R TOA (560) ≥ 0.4 does not guarantee discarding pixels containing sub-pixel clouds.In the ART and FAR codes the second criterion ξ ≤ 1.16 is included to solve this problem at least partially.The estimations show that the use of this criterion leads to discarding pixels where clouds with the reflectance more than 0.2 cover 20% and more of the pixel area.
In our algorithms additional precaution against a bias of the retrieved AOT value due to the influence of clouds or of bright adjacent pixels is provided by the following procedure.The box of neighbours "good" pixels of the same type ("land" or "water") is set up.The box size depends on the spatial resolution of the deployed satellite instrument, number of pixels in the processing image, the scenario and the problem under consideration.For this box the average value of RTOA (λ) is calculated.In doing so 20% of the pixels with minimal values of R TOA (665 nm) and 30% of the pixels with maximal values of R TOA (665 nm) are discarded from the data to be averaged.The determined average value of RTOA (λ) is ascribed to the central pixel in the box.
This averaging procedure implies that the optical characteristics of the aerosol layers are practically the same for all pixels in the compiled pixel box, i.e. the scale of the spatial changes of the aerosol parameters is much larger than the pixel size.It decreases the effects of the random errors of the measurements and of the small-scale variations of the surface albedo.Displacement of the pixel box by one column or one row in the array of the image frame provides the retrieval of the moving-average value of the AOT in the retrieval procedure.Discarding pixels with maximal and minimal values of R TOA (665 nm) in the compiled pixel box from the averaging procedure allows one to eliminate or at least to decrease the effect of pixels which either partially contain clouds (maximal values of R TOA (665 nm)) or include cloud shadows (minimal values of R TOA (665 nm)).
Thus, the pixels with light soil (sand), snow and cloud are discarded, i.e. they are not included in the estimation of the average value of RTOA (λ) for the box.Nevertheless, these pixels can affect signals registered from adjusted pixels because of spreading of the scattered light in the atmosphere.The problem of the influence of adjacent pixels is well-known.Estimations given in (Dave, 1980;Kaufman, 1985) has shown that the influence of adjacent pixels on the signal registered from the chosen pixel extends at the distance of about 5 km because of molecular and aerosol scattering in the atmosphere.Nowadays this feature is taken into account in the procedures of the in-flight calibration of satellite sensors, and particularly for the in-flight calibration of sensors with a high spatial resolution (Richter, 1997).Let us underline that for the AOT retrieval from satellite data with the spatial resolution of about 1 km the small-scale variations of the surface albedo with small spatial fluctuations are efficiently averaged and do not make a noticeable influence.The large-scale spottiness with high albedo (clouds, snow) may affect a measured values of the AOT.In our ART and FAR algorithms (like in the MODIS algorithm (Kaufman et al., 1997) this undesirable effect is somewhat softened by discarding 30% of the pixels with maximal values of R TOA (665 nm) from the data to be averaged over a pixel box.But this procedure may occur insufficient.The AOT values retrieved with the ART using MERIS data for regions of Belsk (Poland) (Fig. 24, left part) and Zvenigorod (Russia) (Fig. 25, left part) averaged over a circle with the radius of 20 km with the centre in the pixel where AERONET measurements were performed are depicted at the left figures.Different signs are given for different cloud situations.We distinguish the situations with numbers of pixels identified as cloudy less than 10%, in a range of 10-50%, and more than 50% in the considered circle with the radius 20 km.The lines show the linear regression between AERONET data and the retrieved AOT values for the corresponding cloudy situations.Table 1 shows the coefficients of the linear regression.As seen, the more nearby the cloudy pixels are, the higher the retrieved AOT values and the AOT dispersion.
The most likely cause of this tendency is the influence of the cloudy pixels (to be more explicit, the pixels identified as cloudy) onto signals from adjusted pixels through the scattering in the atmosphere.To check this assumption let us estimate the scale of "the effect of adjacent cloudy pixels".
This scale is determined by the dispersion of the Point Spread Function (PSF) for the light propagating through the atmosphere.Just this value determines the effective square of land area that contributes to the signal registered from a pixel by a satellite instrument.The area where the contribution from adjacent pixels is noticeable is practically the same for the case of the contribution of adjacent cloudy pixels or from land pixels (see Appendix A).
The dispersion of the PSF for the beam propagating through atmosphere due to aerosol and molecular scattering is estimated in Appendix A. As it follows from estimations the presence of a cloud (snow, sand with high albedo) in some pixel affects signals registered from the nearby pixels up to the distance of about 1 km because of the scattering in aerosol atmosphere and at a few kilometres because of molecular scattering.This influence is expected to be more pronounced in the shortwave part of the spectrum with the increasing contribution of the molecular scattering.
With the goal to decrease the bias of the retrieved AOT due to the impact of the nearby cloudy pixels in the signal from the pixel under retrieval, we included discarding not only cloudy pixels but nearby pixels as well in the retrieval algorithm.The obtained data are given in right parts of Figs.24 and 25.The coefficients of the linear regression between Let us note that this error is not less than, for instance, the errors due to the use of approximate equations in the proposed express-algorithm FAR.
As it is seen from Figs. 24 and 25, the proposed method to decrease the bias of the retrieved AOT values due to effect of nearby cloudy pixels by discarding pixels adjacent to this cloudy pixel, is efficient enough.But it does not eliminate this effect completely.Beside this technique leads to some decrease of the number of pixels for retrieval.The further development of the way to eliminate this bias is desirable.

Conclusions
We have presented in this paper a new FAR technique for the acquisition of the atmosphere spectral AOT and spectral reflectance of the underlying surface from MERIS data.This FAR code may be considered as a faster version of the previously published ART technique (Katsev et al., 2009).Both these codes include radiative transfer computations in the process of the retrieval instead of the LUT and use the least-squares method to determine the AOT and Angström exponent.The basis of both approaches is the accurate and extremely fast radiative transfer code RAY for simulation of the radiative transfer in the atmosphere-underlying surface system with an account for the light polarization.The FAR code keeps all advantages of the ART technique.The main difference between ART and FAR techniques is in the radiative transfer procedures.To speed up satellite data processing and to provide an operational mode, the FAR code uses the combination of RAY calculations and analytical solutions of the radiative transfer theory.The FAR spectral AOT retrieval is carefully checked with ART retrievals and with available ARONET data and found to be satisfactory.Beside FAR is shown to be an efficient and accurate tool for the atmosphere correction and retrieval of the land spectral albedo.It is important that the FAR noticeably speeds up data processing in comparison with the ART In this paper, we introduced the MSA that merges the truncation of the phase function with simple solutions given in Sobolev (1975).The MSA accuracy is checked carefully.This approximation is used only for the computations of the radiative characteristics of the boundary layer of troposphere.
Two more problems of the AOT retrieval are considered.First, we proceeded with estimations of the influence of a priori information about aerosol in the boundary layer on the results of the AOT retrieval from spectral satellite data.The second problem is the influence of clouds in adjacent pixels on the retrieved AOT from cloudless pixels.As for our knowledge, this problem has not been given required attention up to now.The bias in the retrieved AOT due to the clouds in adjacent pixels is confirmed using satellite data.With the goal to decrease this bias, we included discarding not only cloudy pixels but pixels adjacent to the cloudy pixels as well.This procedure decreases the AOT bias noticeably.Nevertheless, the further development in this direction is recommended.
where ω aer is the single scattering albedo for aerosol.
For the molecular atmosphere it can be taken θ 2,mol ≈ 1 while for the aerosol atmosphere with Continental aerosol θ 2,aer ≈ 0.7.
Let us consider atmosphere with τ aer < 1 and τ mol < 1.Then we have from Eqs. (A6) and (A8): i.e. the PSF dispersion is practically independent of the atmospheric optical thickness.It means that the area with adjacent pixels that affects the signal from the pixel under observation is practically the same for clouds and underlying land surface.
In the case of the satellite sensor, when α mol H 1,α aer H 1, from (A9) and (A11), we get: Figure A1 shows dispersions V mol (H ) and V aer (H ) as functions on a sensor altitude H .As seen, in the case of observations from a helicopter (H = 500 m) the effective area of an underling surface that contributes in the registered signal because of the scattering in the aerosol atmosphere is about 0.03 km 2 , while such an area is about 0.3-0.5 km 2 at observation from aircraft at the flight altitude (H = 2−3 km) and achieves 1 km 2 for observation from satellite or aircraft at high altitudes.The molecular scattering in atmosphere leads to increasing the effective surface area that contributes to the registered signal from considered pixel and this influence is expected to be more pronounced in the shortwave part of the spectrum by increasing the contribution of the molecular scattering.
Thus, the presence of a cloud (snow, sand with high albedo) in some pixel affects signals registered from the nearby pixels up the distance of about 1 km because of the scattering in aerosol atmosphere and at a few kilometres because of molecular scattering.Note, that fulfilled estimations are in a good agreement with results reported in (Dave, 1980;Kaufman, 1985) and with requirements to the size of an area for the calibration of satellite optical sensors.Let us emphasize that the PSF dispersion is determined only by the scattering in atmosphere and, as a direct consequence, practically does not depend on the satellite altitude.

Fig. 4 .
Fig.4.Errors of the MSA computations of the reflectance from layers with Continental aerosol at normal incidence dependent on the view zenith angle at the layer optical thickness equal to 0.2 (top), 0.5 (middle), 1.0 (bottom).Data for different cut angles are given in red (0 • ), green (45 • ) and blue (66 • ).

Fig. 12 .
Fig. 12. Correlation between values of AOT at 412.5 nm retrieved with the FAR and ART codes.

Fig. 13 .
Fig. 13.The same as in Fig. 12, except for the wavelength of 550 nm.

Fig. 14 .Fig. 15 .
Fig. 14.Histograms of the distributions of the AOT at the wavelength 412.5 nm retrieved with the FAR and ART codes from MERIS on ENVISAT data for the areas marked in Fig. 10.

Fig. 16 .
Fig. 16.AOT retrieval at 440 nm with codes FAR (top) and ART (center) in comparison with the measurements of AERONET stations in Zvenigorod (Russia, near Moscow) and Belsk (Poland) for periods from March to September in 2008 and 2009 and correlation of values of AOT retrieved with ART and FAR technique (bottom).

Fig. 20 .Fig. 21 .
Fig. 20.Comparison of the AOTs (550 nm) retrieved by the ART algorithm using Continental and experimental phase functions.r is the correlation coefficient.

Fig. 22 .
Fig.22.Spectral reflectance at the TOA computed using code RAY with the retrieved AOT values and atmosphere models used for the retrieval for the AOT equal to 0.3 (left) and 1.0 (right) at λ = 550 nm.

Fig. 23 .
Fig. 23.Surface albedo retrieved with the experimental phase function versus the same value retrieved with the Continental phase function at wavelengths of 412.5 nm (left) and 560 nm (right).r is the correlation coefficient.

Fig. 24 .Fig. 25 .
Fig. 24.The satellite-derived AOT values at λ = 440 nm for the region of Belsk (Poland) versus AERONET AOT.Different signs are given for cloud situations with numbers of pixels identified as cloudy less than 10%, in the range of 10-50%, and more than 50% (see the legends).The lines show the linear regressions between the retrieved AOT values and AERONET data.MERIS data are averaged over a circle with the radius of 20 km with the centre in the pixel where AERONET measurements were performed (left picture).The same but with the additional discarding of pixels in the square 5 × 5 pixels with a cloudy pixel in the centre (right picture).
Fig. A1.Dependence of the dispersions V mol (H) and V aer (H) on the sensor altitude.
The reflectance for the aerosol layer computed with RAY at different truncation angles given in the legend.Aerosol model Continental, the normal incidence; τ = 0.2 (left) and 0.5 (right).

Table 1 .
Coefficients of the linear regression between AERONET data and the retrieved AOT values for different situations.

Table 2 .
Coefficients of the linear regression between the AERONET data and the retrieved AOT values with additional discrimination of pixels adjacent to the cloudy pixel.

Table 2 .
Thus, the only difference between the left and right parts of Figs.24 and 25 is that values of AOT presented in right parts are obtained with the additional discarding of the pixels in the square 5×5 with the centre in the cloudy pixel.Of course, this procedure decreases the number of used pixels in a box.As it was expected, the AOT values retrieved with this additional discarding decrease in comparison with data obtained without this additional discrimination and show better agreement with AERONET data.This correction is more pronounced for the situations with larger cloud fractions.Thus, one can see that the presence of clouds in the nearby pixels could provide noticeable bias of the retrieved AOT values, namely provide overshot AOT value.