Acceleration of radiative transfer model calculations for the retrieval of trace gases under cloudy conditions

Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Acceleration of radiative transfer model calculations for the retrieval of trace gases under cloudy conditions Dmitry S. Efremenko a,n, Diego G. Loyola a, Robert J.D. Spurr b, Adrian Doicu a a b

German Aerospace Centre (DLR), Remote Sensing Technology Institute (IMF), Oberpfaffenhofen, 82234 Wessling, Germany RT Solutions Incorporated, Cambridge, MA, USA

a r t i c l e in f o

abstract

Article history: Received 27 September 2013 Received in revised form 19 November 2013 Accepted 20 November 2013 Available online 27 November 2013

In the independent pixel approximation (IPA), radiative transfer computations involving cloudy scenes require two separate calls to the radiative transfer model (RTM), one call for a clear sky scenario, the other for an atmosphere containing clouds. In this paper, clouds are considered as an optically homogeneous layer. We present two novel methods for RTM performance enhancement with particular application to trace gas retrievals under cloudy conditions. Both methods are based on reusing results from clear-sky RTM calculations to speed up corresponding calculations for the cloud-filled scenario. The first approach is numerically exact, and has been applied to the discrete-ordinate with matrix exponential (DOME) RTM. Results from the original clear sky computation can be saved in the memory and reused for the non-cloudy layers in the second computation. In addition, for the whole-atmosphere boundary-value approach to the determination of the intensity field, we can exploit a ’telescoping technique’ to reduce the dimensionality (and hence the computational effort for the solution) of the boundary value problem in the absence of Rayleigh scattering contributions for higher azimuthal components of the radiation field. The second approach is (for the cloudy scenario) to generate a spectral correction applied to the radiation field from a fast two-stream RTM. This correction is based on the use of principal-component analysis (PCA) applied to a given window of spectral optical property data, in order to exploit redundancy in the data and confine the number of fullstream multiple scatter computations to the first few EOFs (Empirical Orthogonal Functions) arising from the PCA. This method has been applied to the LIDORT RTM; although the method involves some approximation, it provides accuracy better than 0.2%, and a speed-up factor of approximately 2 compared with two calls of RTM. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Fast radiative transfer model Cloud Ozone retrieval

1. Introduction New-generation atmospheric composition sensors such as the Sentinel 5 Precursor, Sentinel 4 and Sentinel 5 [1] will deliver huge amounts of data, with an expected increase of two orders of magnitude in data flow compared to that from

n

Corresponding author. Tel.: þ49 176 39504981. E-mail addresses: [email protected], [email protected] (D.S. Efremenko). 0022-4073/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jqsrt.2013.11.014

current missions. Fast and accurate radiative transfer models (RTMs) are required for systematic processing of forthcoming Sentinel measurements. Operational systems required to process these data will need to be optimized in order to cope with near-real-time data turn-over requirements. The presence of clouds is an important aspect of trace gas retrievals from such sensors. Ozone retrieval from the GOME-2 instrument, for example, is performed in two stages [2]. In the first step, ancillary cloud parameters (such as cloud fraction c, cloud top height, and cloud optical thickness) are obtained using the OCRA (Optical

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

Cloud Recognition Algorithm) and ROCINN (Retrieval of Cloud Information using Neural Networks) algorithms [3]. Then, in the main algorithm for retrieval of ozone column densities from UV backscatter measurements in the 325– 335 nm window, this cloud information is a necessary ingredient in the forward model RTM simulations for satellite intensities and associated weighting functions [4]. 2-D or 3-D RTMs are far too time-consuming for operational data processing. Existing retrieval techniques are entirely one-dimensional, assuming horizontal homogeneity on the pixel scale. In the independent-pixel approximation (IPA), radiation fields are computed on a pixel-by-pixel basis with plane-parallel scattering RTMs. In the IPA, cloud fraction c is the percentage of each pixel in satellite imagery that is covered with clouds, and the radiance for the current pixel is computed as a superposition of solutions for a clear sky case and a cloudy case: I ¼ cIcloud þ ð1  cÞI clear

ð1Þ

where I cloud is the radiance for a fully cloudy pixel and I clear is the radiance for a cloud free pixel. Since the cloud scenario consists of one to three contiguous cloud layers embedded in an otherwise clear-sky atmosphere, calculations for I cloud and I clear will differ only for these cloud layers, and it makes sense to save preliminary RTM results from the clear-sky computation and reuse them in the cloudy-scenario computation. Recently, optical property dimensionality reduction techniques for accelerated radiative transfer performance have been been applied to remote sensing total ozone retrievals [5–7]. The dimensionality of the hyper-spectral input optical property data set is reduced by means of principal component analysis (PCA) [8,9] or other analogous locally linear transformations which exploit redundancy in the optical data. Full multi-stream RTM computations are then calculated for a reduced input optical data space, and the results are used to develop a correction which is then applied to approximate intensity output from a fast twostream model. For more information on the technique, see [10]. In this paper, we use a two-stream model to compute a correction factor which can be applied to the cloudy computation in order to avoid an additional call to the full-stream RTM. This paper is structured as follows. In Section 2, we briefly summarize the discrete ordinate method for solving the radiative transfer equation (RTE) in the presence of cloud layers. In Section 3, we describe the acceleration technique based on storage and reuse of clear-sky computations, plus the use of the boundary-value problem ‘telescoping technique”, as applied to the DOME matrix-operator model. Section 4 contains a summary of performance enhancement using Principal Component methods, and is mainly concerned with the development of cloud-intensity corrections based on the use of fast two-stream model calculations as applied to clear sky computations using the LIDORT RTM. The paper concludes with a short summary.

59

here will put our explanations in a proper context. The RTE for diffuse radiance in the absence of polarization is Z dI ωðτÞ μ ðτ; ΩÞ ¼  I ðτ; ΩÞ þ P ðτ; Ω; Ω′ÞI ðτ; Ω′Þ dΩ′ dτ 4π 4π ωðτÞ F 0 e  τ0 ðτÞ P ðτ; Ω; Ω0 Þ ð2Þ þ 4π Here τ is the optical depth of the medium, ωðτÞ the single scattering albedo, F0 the incident solar flux, P the scattering phase function, τ0 the solar optical depth, and qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi Ω0 ¼ ð 1 μ20 cos φ0 ; 1  μ20 sin φ0 ;  μ0 Þ with fμ0 ; φ0 g the incident solar direction, where μ0 4 0 in this definition. The quantity Ω is defined similarly in terms of the polar stream cosine μ and azimuth angle φ. Our convention is that μ 40 for upward directions and μ o0 for downward directions. In the pseudo-spherical model, attenuation of the direct solar beam is treated for a spherically curved atmosphere, while the scattering is treated in a planeparallel geometry [4]. Eq. (2) is discretized in the angular domain by first considering an azimuthal expansion of the radiance field: M

Iðτ; ΩÞ ¼ ∑ I m ðτ; μÞ cos mðφ  φ0 Þ;

ð3Þ

m¼0

and second by defining a set of N do Gaussian quadrature points and weights fμk ; wk gk ¼ 1;Ndo for the computation of the multiple-scattering integral in the polar direction. Here Ndo is the number of discrete ordinates in the polar hemisphere, and M is the number of azimuthal harmonics. For an inhomogeneous atmosphere, we consider a spatial discretization of the atmosphere with N levels, i.e., τ1 4 τ2 4 ⋯ 4 τN (Fig. 1). A layer j is bounded above by the level τj and below by the level τj þ 1; the number of layers is N  1. The single scattering albedo and the phase function are assumed to be constant within each layer, and for the layer j the optical thickness is Δτj ¼ τj τj þ 1 . The above discretization process leads to a set of coupled linear differential equations for the radiance field in the layer j: dim ðτÞ ¼ Amj im ðτÞ þ bmj ; dτ þ



τj þ 1 r τ rτj ;

ð4Þ

7

where im ðτÞ ¼ ½im ðτÞ im ðτÞT with ½im ðτÞk ¼ I m ðτ; 7 μk Þ is the radiance vector in discrete ordinate space, Amj is a layer matrix which encapsulates the scattering properties of the

2. RTM under cloudy conditions Although the discrete ordinate solution of the 1-D RTE is well known [11,12], a short mathematical exposition

Fig. 1. The atmosphere, discretized in N  1 homogeneous layers.

60

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

Table 1 Organization of the whole-atmosphere global matrix Um and source vector Bm in Eq. (8). Row\column

2 N do

2 N do

2 N do

…

2 N do

2 N do

½0; I

0 U2m;1

… …

0 0

0 0

0 Bm;1

2N do

U1m;1 0

0 0 U2m;2 … 0

0

0

Bm;2

… 0

U1m;2 … 0

…

… 2N do

… …

…

…

0

0

0

…

U2m;N  1 ½I; Rm 

… Bm;N  1

N do

U1m;N  1 0

N do 2N do

j1

Bmj  exp  ∑ Δτk :

ð9Þ

k¼1

Bm

Um

thicknesses of layers above j, that is !

2. The matrix operator method is based on the interaction principle [15]. Eq. (5) is rewritten as " þ # " þ þ #" i  # " þ # Σ mj imj R mj Tmj mj ¼ þ ð10Þ    þ  ; Σ mj T R imj þ 1 imj þ 1 mj mj

rm

layer, and bmj is a layer vector. The solution of the “layer” equation (4) involves the following steps: 1. derive the integral form of equation (4), which relates the radiances at the boundary levels imj and imj þ 1 , 2. perform a spectral decomposition of the layer matrix 1 Amj , i.e., Amj ¼ Vmj Λmj Vmj , where Vmj is the eigenvector matrix and Λmj is the eigenvalue matrix, 3. integrate over the layer vector bmj by employing the average secant approximation for the solar beam attenuation [4], and finally, 4. apply a scaling transformation [13] to avoid the numerical instability.

þ þ and Tmj are the reflection and the transmiswhere R mj 7 sion matrix, respectively, and Σ mj is the source vector [16]. The computation process is organized recursively by using the concept of a “stack”. Starting from the top of the atmosphere, the interaction principle and the adding algorithm are used to unify the first and the second layer into a single layer. The resulting layer, with effective reflection and transmission matrices and effective source function, is called “a stack”. The procedure is then repeated until the last layer is added to the stack. Finally, we are led to a matrix equation for the stack of all atmospheric layers as in (10).

3. Performance enhancement for DOME in the IPA 3.1. General considerations

The final expression for algebraic layer equation takes the form: U1mj imj þU2mj imj þ 1 ¼ Bmj ;

ð5Þ

while the boundary conditions at the top and bottom of the atmosphere are ½0; Iim1 ¼ 0;

ð6Þ

and ½I; Rm imN ¼ rm ;

ð7Þ

respectively. In (6) and (7), I is a identity matrix, Rm is the reflection matrix, and rm is the reflection vector [14]. Solving for the radiation field in a multi-layer medium can be done via the whole-atmosphere approach, or via the matrix operator method:

1. In the whole-atmosphere approach, layer equations (5) are assembled into a global linear matrix algebra problem for the entire atmosphere: Um im ¼ Bm ;

ð8Þ

where layouts of the global matrix Um and source vector Bm are shown in Table 1. The global matrix has a banded structure, and it can be stored in compressed form and inverted using standard methods. It is important to note that matrices U1mj and U2mj do not depend on the optical parameters for layers other than j, but source vector Bmj depends exponentially on optical

For the retrieval of atmospheric constituents from satellite-measurements of Earth-shine radiation, we have developed a software tool at the German Aerospace Centre (DLR), which includes the DOME RTM based on the discrete ordinate method with matrix exponential formalism [14]. DOME solves for the complete radiance field using the whole-atmosphere approach or the matrix operator method. For computing the radiance field at specific directions away from the discrete ordinates, DOME uses either the false discrete ordinate method [17–19] or the source function integration technique [20,21]. In this paper, a clear-sky scenario is defined to be a Rayleigh-scattering atmosphere with molecular absorption but with no aerosols. In such a scattering regime, clear sky computations are performed only for the m ¼0, 1, 2 azimuthal harmonics, since the azimuthal expansion coefficients of the (Rayleigh) phase function vanish for m 42. Phase functions for cloud scattering are typically strongly peaked in the forward direction, and will require many more azimuthal harmonics to capture the scattering accurately. The key idea for optimizing the RTM computation for an atmosphere containing a cloud layer is to note that the inclusion of a cloud layer into a clear-sky atmosphere does not change the optical properties of the clearsky layers. Thus, it is helpful to store any solution results derived from the call to the RTM for the clear sky scenario, and use these results again for the call to the scenario with cloud layer(s) present. In addition, for a cloudy scenario with one or more contiguous cloud layers in an otherwise clear-sky atmosphere, it is helpful to use the so-called telescoping technique [22] to speed up the computation. Since there

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

is no scattering for azimuthal harmonics m 42 in the ‘Rayleigh’ layers, matrices Amj are diagonal in these layers, and the corresponding vectors bmj are absent. In this case, it is only necessary to solve the ’global’ boundary value problem for the few layers containing clouds, and compute radiances at the remaining levels recursively through the Beer–Lambert–Bouguer. Analysis of the speedup due to the telescoping technique can be found in [23]. 3.2. Optimization for the whole-atmosphere approach The insertion of a cloud layer in a clear sky atmosphere will affect the atmospheric layering scheme. This depends on the cloud-top height and the cloud geometrical thickness, and the possible options are illustrated in Fig. 2. The simplest case involves a cloud with optical thickness Δτ introduced into the layer j0 which has clear-sky optical depth Δτ0 (Fig. 2b). In this case, when solving the clear-sky problem we store the matrices U1mj and U2mj for all layers j aj0 , as well as the source vectors Bmj for all layers j aj0 . When solving the cloudy-sky problem, we use the clear-sky layer equation (5) for all layers j o j0 ; for the layers j 4 j0 we take account of Eq. (9) and use a source vector corrected for the change in attenuation: Δτ0  Δτ : Bcor mj ¼ Bmj e

ð11Þ

A similar correction is applied to the reflection vector: Δτ0  Δτ : rcor m ¼ rm e

ð12Þ

If the boundary of the cloud splits a layer, as shown in Fig. 2c, then these corrections are applied to clear-sky layers situated above and below the split layers. For the split layers themselves, we must store the solutions to the homogeneous RTE obtained for the original clear-sky calculation.

61

 as in (11), multiply the source vectors Σ m7 of stack 2 by expðΔτ0  ΔτÞ.

Thus, for a cloudy atmosphere, we must derive the matrix equation (10) for the cloud layer itself, and then apply the interaction principle and the adding algorithm for the three layers depicted in Fig. 3b. If the cloud boundary splits the layer j, then the þ þ reflection and the transmission matrix R mj and Tmj , respectively, will be changed. Unfortunately, it is not possible to deploy a simple analytical expression for the correction of these matrices; the best we can do is store RTE solution vectors for the split layers, and from these reconstruct the two matrices from scratch. 3.4. Performance enhancement for DOME To establish the efficiency of the above optimization in the context of the total ozone retrieval problem for the spectral domain 325–335 nm with 88 spectral points, we discretize the atmosphere with height resolutions of 1 km between 0 and 25 km, 2.5 km between 25 and 50 km and 10 km between 50 and 100 km. The total number of layers is 40. Calculations are performed for a standard atmosphere [24] with an ozone profile of total content 300 DU taken from the column classified TOMS Version 8 climatology [25]. Rayleigh scattering properties are taken from the parametrization in [26]. The (sun-normalized) radiance at the top of the atmosphere is computed for a viewing zenith angle (VZA) 601 and a nadir viewing geometry (VZA¼0), a solar zenith angle of 301, and surface albedo of 0.1.

3.3. Optimization for the matrix operator method In the matrix operator method, individual layers are assembled into one multi-layer stack with effective scattering properties. Assuming that the layer j0 contains a cloud, the optimization process can be organized as follows:

 store in stack 1 all upper layers j ¼ 1; …; j0  1,  store in stack 2 all lower layer j ¼ j0 þ 1; …; N,

Fig. 3. Position of a cloud between two stacks.

Fig. 2. Cloud position with respect to atmosphere layers. (a) a "clear" atmosphere, (b) a cloud occupies the whole layer and (c) a cloud splits the layers.

62

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

Table 2 Computational times for DOME IPA computations with and without optimization. MOM refers to the matrix operator method, while WA refers to the whole atmosphere approach. VZA ɛ (deg)

0 0 60 60

M MOM (s)

10  4 5 10  6 8 10  4 6 10  6 10

24 26 25 27

MOM opt. (s)

Relative speedup (%)

WA (s)

WA opt. (s)

Relative speedup (%)

14 15 15 15

92 86 83 90

21 22 21 24

11 12 12 13

95 91 87 92

We consider a cumulus cloud with a modified-gamma particle size distribution:   γ  α a pðaÞ p aα exp  ; γ amod and a droplet size ranging between 0.02 and 50.0 μm [27]. The cloud top height and the cloud geometrical thickness are htop ¼ 5 km and Δh ¼ 1 km, respectively, while the cloud optical thickness is Δτ ¼ 1. In our simulations we use N do ¼ 16 discrete ordinates in the polar hemisphere. The solution in the VZA directions is obtained by using source function integration. To achieve convergence, we compute azimuthal harmonics including up to M¼ 5 for an accuracy of ɛ ¼10  4, and up to M¼8 harmonics for accuracy of ɛ ¼10  6 for VZA¼01 and M¼6 for an accuracy of ɛ ¼ 10  4, and up to M ¼10 harmonics for accuracy of ɛ ¼10  6 for VZA¼601. In the cloudy scene computation, solutions for azimuthal modes m 4 2 are computed by using the telescoping technique. Computational times for these cases, as well as the relative speedup, are shown in Table 2. The computations are performed on an Intel Xeon E5-1620 3.60 GHz machine, with OpenSUSE 12.2 Linux operating system. With the telescoping technique, computations for m 42 are relatively fast. For the m¼ 0, 1, 2 harmonics, the only additional computations for the cloudy-sky problem are the spectral decomposition (eigensolution) of the matrix Amj for the cloud layer and the subsequent adjustments needed to solve for the complete field. A combination of these considerations yields nearly double the speed-up.

4. LIDORT performance enhancement in the IPA using PCA methods 4.1. PCA-performance enhancement for LIDORT LIDORT is a multiple-scattering scalar RTM that uses the multi-stream discrete ordinate method for solving the RTE plus a whole-atmosphere boundary-value problem approach to determine the complete radiance field in a multi-layer optically stratified atmosphere [28]. LIDORT uses the pseudo-spherical approximation, treating solarbeam attenuation for a spherically curved atmosphere (scattering is plane-parallel). Although it is possible to run LIDORT in two-stream model (with one stream in the polar hemisphere), the performance enhancement to be described

below was done using a stand-alone multiple-scatter twostream code [29] that is fully compatible with LIDORT itself. In order to enhance performance, it is helpful to greatly reduce the number of full multiple scatter calculations using a combination of a fast two-stream RTM and a much smaller set of optical properties. In this regard, one can apply principal component analysis (PCA) to the meanremoved correlation matrix associated with the space of optical properties to be used for a given window. If the optical data are strongly correlated, the first few principal components will capture the vast majority of the data variance, and it is then only necessary to carry out full calculations for a small number of optical profiles. In practice, one carries out full LIDORT and fast two-stream RT computations for these first few optical-property EOFs (Empirical Orthogonal Functions), and the differences between these results are then used to develop a wavelength-dependent correction factor which is applied to fast two-stream computations carried out for the original full-spectrum optical data set. This was the method used in [30]. In [6], single scatter computations were done separately, with correction factors limited to the multiple-scatter part of the field. Also in [6,7], the PCA-based performance enhancement method was given a linearization facility for fast generation of weighting functions as well as intensity; this method was then applied to total ozone retrieval in the 325–335 nm fitting window. It was found that for 88 spectral points, it is possible to obtain ozone retrieval accuracy better than 0.5% with just 2 EOFs, with a speed-up of 3.5 for eightstream (full sphere) LIDORT calculations. In [7], a number of other dimensionality reduction techniques were applied to linearized RT calculations for the above total ozone problem. We also note the VLIDORT PCA-based application by [31], which used a second-order-of-scattering (2OS) model to develop a fast correction for the polarization signature in the O2 Aband. In [6], the PCA performance enhancement is applied to the scalar intensity, and in the total ozone retrieval [32] a polarization correction factor is applied to this intensity – this factor is supplied in the operational algorithm through use of a look-up table [32,33]. Since polarization effects are important in such a Rayleigh-atmosphere situation, a related investigation is currently underway to use the PCA method to enhance performance for the full Stokes intensity, without the need for an additional LUT. 4.2. Two-stream correction for the cloud computation As noted above, the main idea is to compute a twostream correction factor f ðλÞ based on two-stream and multi-stream model calculations for the reduced space of optical parameters obtained from PCA.The multi-stream solution IðλÞ is then approximated as IðλÞ  I TS ðλÞf ðλÞ; TS

ð13Þ

where I ðλÞ is the two-stream solution. In the following we will develop a suitable correction factor to speed-up the radiative transfer calculation in a cloudy atmosphere. First, we propose the following computational formula for the multi-stream solution in a

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

cloudy atmosphere: I cloud ðλÞ 

I TS ðλÞ K ðλÞ: I clear ðλÞ cloud I TS clear ðλÞ

ð14Þ

Here I TS cloud ðλÞ is the two-stream solution for the cloudy scenario, I TS clear ðλÞ is the two-stream solution for clear-sky scenario and KðλÞ is the correction factor to be determined. Second, applying the dimensionality reduction techniques for computing the multi-stream solution for a clear sky, i. e., I clear ðλÞ  I TS clear ðλÞf ðλÞ, (14) becomes I cloud ðλÞ  I TS clear f ðλÞ

I TS cloud ðλÞ I TS clear ðλÞ

K ðλÞ

ð15Þ

With this approach, the computation of the multi-stream solution for a cloudy sky I cloud requires only one additional call to of the two-stream solution I TS cloud for each spectral point. 4.3. Approximations for the correction factor and numerical simulations Correction factors KðλÞ are pre-computed for various values of the cloud parameters using the following inverted form of Eq. (15): K ðλÞ ¼

I cloud ðλÞ I TS clear ðλÞ : I clear ðλÞ I TS cloudy ðλÞ

ð16Þ

Next, KðλÞ is linearly interpolated in the spectral domain as KðλÞ  kðλ  λÞ þ b

The 8 test cases are outlined in Table 3. The first 5 cloud cases are based on particle size distribution parameters taken from [27], with Mie calculations for the phase function. Test cases 6–8 are based on the use of the Henyey–Greenstien phase function [34] for asymmetry parameter values as indicated. In all cases, the cloud single scattering albedo was set to 0.99. Apart from these cloud specifications, the clear-sky optical data is the same as in the previous section on DOME optimization, with the sole exception of the ozone profile which was set to a total column of 450 DU. For test case 1, the ratios of the radiances computed with two- and multi-stream models are shown in Fig. 4 while the corresponding correction factors KðλÞ are plotted in Fig. 5. Results corresponding to test case 5 are illustrated in Figs. 6 and 7. Since the two-stream and the multi-stream solution have broadly similar spectral dependence, the correction factor KðλÞ is close to linear with wavelength, although strictly speaking KðλÞ is a slightly non-monotonic function of λ. For test cases 1–5, the linear spectral correction (14) yields a maximum error of about 0.1–0.2 %, which is well within accuracy requirements for the total ozone retrieval.

1.1

Full-stream / two-stream, cloudy Full-stream / two-stream, clear sky 1.05

ð17Þ

where k and b are constants, the values of which are stored in look-up tables, and λ is the mean wavelength for the spectral window (in this case λ ¼ 330 nm). Substituting (17) in (14) gives

1

0.95



I TS ðλÞ  k λλ þb : I cloud ðλÞ  I clear ðλÞ cloud TS I clear ðλÞ

ð18Þ

To investigate the accuracy of this approximation, we performed simulations for different types of clouds and different values of the cloud geometrical thickness (CGT), cloud top height (CTH), and cloud optical thickness (COT). Apart from the clouds, the atmospheric layering and optical scenario are the same as that used in Section 3.4.

0.9 325

Test case

CTH (km)

CGT (km)

COT Cloud type Average Maximum error (%) error (%)

1 2 3 4

4 3 3 3

2 1 1 1

1 1 5 5

5 6

9 20

1 15

20 10

7

30

20

20

40

30

30

Cum. clean Cum. clean Cum. clean Cum. polluted Cum. clean HG g ¼ 0:9 HG g ¼ 0:9 HG g ¼ 0:9

0.01 0.01 0.05 0.2

0.04 0.02 0.2 0.1

0.1 0.5 (0.02) 0.6 (0.03) 0.6 (0.03)

0.2 1.4(0.2)

326

327

328

329

330

331

332

333

334

335

Wavelength [nm] Fig. 4. The ratios of the radiances computed with two- and multi-stream models for the test case 1 from Table 3.

1.0085

Table 3 Input cloud parameters and approximation errors for Eq. (14). Numbers in parenthesis are errors for the non-linear approximation in Eq. (19).

8

63

Correction factor K(λ) Linear approximation

1.008

1.0075

1.007

1.0065

1.4(0.2)

1.006 325

326

327

328

329

330

331

332

333

334

335

Wavelength [nm] 1.6(0.2) Fig. 5. The correction factor KðλÞ for test case 1 from Table 3 and its linear approximation.

64

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

1.1

0.865

Full-stream / two-stream, cloudy Full-stream / two-stream, clear sky

0.86 0.855

Correction factor K(λ) Linear approximation Nonlinear approximation

1.05 0.85 0.845 0.84

1

0.835 0.83 0.95 0.825 0.82 0.9 325

326

327

328

329

330

331

332

333

334

335

0.815 325

326

327

0.968

329

330

331

332

333

334

335

Fig. 9. The same as in Fig. 8, but for test case 8.

Fig. 6. The same as in Fig. 4, but for test case 5.

0.97

328

Wavelength [nm]

Wavelength [nm]

about 0.5%. To improve the approximation of KðλÞ and to take into account the quasi-periodic structure of KðλÞ, we introduce an additional term in (17), that is

Correction factor K(λ) Linear approximation

0.966

KðλÞ  kðλ  λÞ þ b þvsO3 ðλÞ:

0.964

where sO3 is the O3 absorption cross section at temperature 270 K, taken from [35] and convolved (with the GOME slit function) to the 88-point wavelength grid in our 325– 335 nm window, and v is a constant. The average and the maximum errors for this non-linear approximation are also indicated in Table 3. This non-linear correction reduces the maximum error to less than 0.2%, while the average error is below 0.03%. Cross sections taken at other temperatures yield very similar levels of accuracy. Figs. 8 and 9 illustrate the results corresponding to test cases 6 and 8, respectively.

0.962 0.96 0.958 0.956 325

326

327

328

329

330

331

332

333

334

335

Wavelength [nm] Fig. 7. The same as in Fig. 5, but for the test case 5.

ð19Þ

5. Conclusions 0.95 0.945 0.94

Correction factor K(λ) Linear approximation Nonlinear approximation

In this paper, we proposed two novel techniques to accelerate the RTM for computations based on cloud-layer scenarios. In particular:

0.935 0.93 0.925 0.92 0.915 0.91 0.905 0.9 325

326

327

328

329

330

331

332

333

334

335

Wavelength [nm] Fig. 8. The correction factor KðλÞ for test case 6 from Table 3, its linear approximation (17) and its nonlinear approximation (19).

In test cases 6–8, which are characterized by large values of the optical and geometrical thicknesses, high cloud tops and high values of the asymmetry parameter, KðλÞ has a marked non-linear dependence. The maximum error of the linear approximation is 1.5%, while the average error is

1. We store and reuse solution quantities obtained from clear-sky simulations. They depend on the spectral decomposition of the layer matrix. For layers below a cloud layer, the source function has to be corrected for the change in attenuation. This approach is implemented in the DOME RTM, for the matrix operator method and the whole-atmosphere approach. The method is exact (no approximations are invoked). 2. In the context of the PCA-based performance enhancement with LIDORT, we approximate the solution for the cloudy atmosphere by a relation computed with a twostream model and a correction factor. The correction factor dependence on wavelength for geometrically thin clouds is almost linear, with maximum radiance error o 0:2% and average error o 0:1%. For geometrically thick clouds, the introduction of a non-linear correction based on the ozone cross section shape yields maximum radiance error of 0.2% and average error below 0.05 %.

D.S. Efremenko et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 58–65

Numerical simulations show that the proposed methods provide almost double speedup of computations. These approaches can also be used for the linearized RTMs, whereby the derivatives of the radiances with respect to specific atmospheric properties are required for retrieval. Linearization of the first approach is straightforward. The second approach requires the derivative of the correction factor KðλÞ. One way to proceed here is to extend the look-up table calculations for KðλÞ for different values of the parameters to be retrieved, and then the derivatives for KðλÞ can be derived by finite difference estimates; however, our computations revealed a weak dependence of KðλÞ on the ozone total column. This topic will be the subject of future papers. References [1] Ingmann P, Veihelmann B, Langen J, Lamarre D, Stark H. Bazalgette Courreges-Lacoste G. Requirements for the GMES Atmosphere Service and ESA's implementation concept: sentinels-4/-5 and -5p. Remote Sens Environ 2012;120:58–69. [2] Van Roozendael M, Loyola D, Spurr R, Balis D, Lambert J-C, Livschitz Y, Valks P, Ruppert T, Kenter P, Fayt C, Zehner C. Ten years of GOME/ERS2 total ozone data: the new GOME Data Processor (GDP) Version 4: I. Algorithm Description. J Geophys Res 2006;111(D14311). [3] Loyola D, Thomas W, Livschitz Y, Ruppert T, Albert P, Hollmann R. Cloud properties derived from GOME/ERS-2 backscatter data for trace gas retrieval. IEEE Trans Geosci Remote Sens 2007;45:2747–58. [4] Spurr RJD. Simultaneous derivation of the intensities and weighting functions in a general speudo-spherical discrete ordinate radiative transfer treatment. J Quant Spectrosc Radiat Transf 2002;75:129–75. [5] Natraj V, Shia RL, Yung YL. On the use of principal component analysis to speed up radiative transfer calculations. J Quant Spectrosc Radiat Transf 2010;111:810–6. [6] Spurr R, Natraj V, Lerot C, Van Roozendael M, Loyola D. Linearization of the Principal Component Analysis method for radiative transfer acceleration: application to retrieval algorithms and sensitivity studies. J Quant Spectrosc Radiat Transf 2013;125:1–17. [7] Efremenko D, Doicu A, Loyola D, Trautmann T. Optical property dimensionality reduction techniques for accelerated radiative transfer performance: Application to remote sensing total ozone retrievals. J Quant Spectrosc Radiat Transf. 2014;133:128–35. [8] Pearson K. On lines and planes of closest fit to systems of points in space. Phil Mag 1901;2:559–72. [9] Jolliffe IT. Principal component analysis. Berlin: Springer-Verlag; 1986. [10] Natraj V, Jiang X, Shia RL, Huang X, Margolis JS, Yung YL. Application of the principal component analysis to high spectral resolution radiative transfer: a case study of the O2A band. J Quant Spectrosc Radiat Transf 2005;95:539–56. [11] Chandrasekhar S. Radiative transfer. New York: Dover Publications; 1960. [12] Stamnes K, Tsay S-C, Wiscombe W, Jayaweera K. Numerically stable algorithm for discrete ordinate method radiative transfer in multiple scattering and emitting layered media. Appl Opt 1988;27:2503–9. [13] Karp AH, Greenstadt J, Fillmore JA. Radiative transfer through an arbitrary thick scattering atmosphere. J Quant Spectrosc Radiat Transf 1980;24:391–406.

65

[14] Doicu A, Trautmann T. Discrete-ordinate method with matrix exponential for a speudo-spherical atmosphere: scalar case. J Quant Spectrosc Radiat Transf 2009;110:146–58. [15] Plass GN, Kattawar GW, Catchings FE. Matrix operator theory of radiative transfer. 1: Rayleigh scattering. Appl Opt 1973;2:314–29. [16] Spurr RJD, Christi MJ. Linearization of the interaction principle: analytic Jacobians in the Radiant model. J Quant Spectrosc Radiat Transf 2007;103:431–46. [17] Evans KF, Stephens GL. A new polarized atmospheric radiative transfer model. J Quant Spectrosc Radiat Transf 1991;46:413–23. [18] Chalhoub E, Garcia R. The equivalence between two techniques of angular interpolation for the discrete-ordinates method. J Quant Spectrosc Radiat Transf 2000;64:517–35. [19] Liu Q, Weng F. Advanced doubling-adding method for radiative transfer in planetary atmosphere. J Atmos Sci. 2006;63:3459–65. [20] Chandrasekhar S. On the radiative equilibrium of a stellar atmosphere. II. Astrophys J 1944;100:76–86. [21] Stamnes K. On the computation of angular distributions of radiation in planetary atmospheres. J Quant Spectrosc Radiat Transf 1982;28: 47–51. [22] Spurr RJD, LIDORT and VLIDORT. Linearized pseudo-spherical scalar and vector discrete ordinate radiative transfer models for use in remote sensing retrieval problems. In Kokhanovsky AA, editor. Light scattering reviews, vol. 3. Berlin: Springer Verlag; 2008. 229–71. [23] Efremenko D, Doicu A, Loyola D, Trautmann T. Acceleration techniques for the discrete ordinate method. J Quant Spectrosc Radiat Transf 2013;114:73–81. [24] U.S. Standard Atmosphere, 1976. U.S. Government Printing Office, Washington, D.C., 1976. [25] McPeters RD, Logan JA, Labow GJ. Ozone climatological profiles for Version 8 TOMS and SBUV retrievals. American Geophysical Union, Fall Meeting 2003, abstract #A21D-0998. [26] Chance KV, Spurr RJD. Ring effect studies: Rayleigh scattering, including molecular parameters for rotational Raman scattering, and the Fraunhofer spectrum. Appl Opt 1997;36:5224–30. [27] Hess M, Koepke P, Schult I. Optical properties of aerosols and clouds. The software package OPAC. Bull Am Meteorol Soc 1998;79:831–44. [28] Spurr RJD, de Haan Johan, van Oss Roeland, Vasilkov A. Discreteordinate radiative transfer in a stratified medium with first-order rotational Raman scattering. J Quant Spectrosc Radiat Transf 2008;109:404–25. [29] Spurr R, Natraj V. A linearized two-stream radiative transfer code for fast approximation of multiple-scatter fields. J Quant Spectrosc Radiat Transf 2011;112:2630–7. [30] Natraj V, Jiang X, Shia R, Huang X, Margolis JS, Yung YL. Application of principal component analysis to high spectral resolution radiative transfer: a case study of the O2A band. J Quant Spectrosc Radiat Transf 2005;95:539–56. [31] Natraj V, Shia R, Yung YL. On the use of principal component analysis to speed up radiative transfer calculations. J Quant Spectrosc Radiat Transf 2010;111:810–6. [32] Roozendael VM, Spurr R, Loyola D, Lerot C, Balis D, Lambert JC, et al. Sixteen years of GOME/ERS2 total ozone data: the new direct-fitting GOME Data Processor (GDP) Version 5: I. algorithm description. J Geophys Res 2012;117(D03305):1–18. [33] Lerot, et al. The GODFIT algorithm: a direct fitting approach to improve the accuracy of total ozone measurements from GOME. Int J Remote Sens 2013;31:543–50. [34] Henyey LG, Greenstein JL. Diffuse radiation in the galaxy. Astrophys J 1941;93:70–83. [35] Daumont D, Brion J, Charbonnier J, Malicet J. Ozone UV spectroscopy. I. Absorption cross-sections at room temperature. J Atmos Chem 1992;15:145–55.