- Email: [email protected]
Adjoint problem of radiative transfer for a pseudo-spherical atmosphere and general viewing geometries
ARTICLE IN PRESS Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
Adjoint problem of radiative transfer for a pseudo-spherical atmosphere and general viewing geometries A. Doicu , T. Trautmann Remote Sensing Technology Institute, German Aerospace Centre, Oberpfaffenhofen, Wessling, Germany
a r t i c l e i n f o
abstract
Article history: Received 14 October 2008 Received in revised form 19 January 2009 Accepted 19 January 2009
We formulate the adjoint radiative transfer for a pseudo-spherical atmosphere and various retrieval scenarios. The single scattering radiance is computed in a spherical atmosphere by using the source integration technique, while for the multiple scattering radiance we formulate an one-dimensional adjoint radiative transfer equation in a plane-parallel atmosphere. The adjoint solution of the radiative transfer equation is obtained by employing the discrete ordinate method with matrix exponential. We provide an abbreviated derivation of our formalism as well as a discussion of the numerical implementation of the theory. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Radiative transfer Pseudo-spherical atmosphere Adjoint operator
1. Introduction The adjoint approach allows the linearization of radiative transfer problem with respect to the optical properties of the atmosphere and the surface. The application of this approach to weighting functions calculation requires the formulation and the solution of the adjoint radiative transfer equation. For a plane-parallel atmosphere, formulations of the adjoint radiative transfer have been given by Marchuk [1], Box et al. [2], Ustinov [3–5], Rozanov and Rozanov [6] and Landgraf et al. [7]. For a pseudo-spherical atmosphere, Walter et al. [8] applied the adjoint approach to a nadir viewing geometry, in which case, the singly scattered radiation is calculated for a spherical geometry, while the multiply scattered radiation is computed for a plane-parallel geometry. Recently, Ustinov [9] extended the general adjoint approach for nadir geometries to limb geometries, but determined the total scattered radiation in a plane-parallel manner. Several solution methods have been employed for solving the adjoint radiative transfer problem. In this context we mention the finite difference method [10], the discrete ordinate method [11,12], and the Gauss–Seidel iteration technique [7]. The goal of our paper is to formulate the adjoint radiative transfer for a pseudo-spherical atmosphere, and for both nadir and limb viewing geometries. The adopted solution method for solving the adjoint radiative transfer problem is the discrete ordinate method with matrix exponential [13]. This method is unconditionally stable for arbitrarily large optical depths and can be simply extended to the present application.
Corresponding author.
E-mail address: [email protected] (A. Doicu). 0022-4073/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2009.01.027
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
465
2. The adjoint radiative transfer approach Depending on the choice of the forward operator and the treatment of boundary conditions, equivalent formulations of the adjoint radiative transfer can be derived [4,6]. In the present analysis we adopt the technique described in [4,7], which assumes the following steps: (1) representation of the boundary-value problem for the radiative transfer in an operator form; (2) derivation of the adjoint radiative transfer operator from the Lagrange identity; (3) derivation of the adjoint radiative transfer problem by using the integral representation of the radiance measured by a satellite instrument. The main results of this section will be written with cursive letters, and for reasons of conciseness, their proofs will be only sketched. 2.1. The forward radiative transfer equation In a pseudo-spherical atmosphere, the boundary-value problem for the diffuse radiance consists in the inhomogeneous differential equation
m
dI ðr; XÞ ¼ sext ðrÞIðr; XÞ þ Jðr; XÞ dr
(1)
the top-of-atmosphere boundary condition (r ¼ r TOA ), Iðr TOA ; X Þ ¼ 0
(2)
and the surface boundary condition (r ¼ r s ), Z A Iðr s ; Xþ Þ ¼ Iðr s ; X Þjm jrðXþ ; X Þ dO þ Q s ðr s ; Xþ Þ.
p
(3)
2p
The source function Jðr; XÞ ¼ J ss ðr; XÞ þ J ms ðr; XÞ includes the contribution of the single and the multiple scattering terms Jss ðr; XÞ ¼ F sun
sscat ðrÞ sun Pðr; X; Xsun Þetext ðjrrTOA jÞ 4p
(4)
and Jms ðr; XÞ ¼
sscat ðrÞ 4p
Z
Pðr; X; X0 ÞIðr; X0 Þ dO0 ,
(5)
4p
respectively, while the surface function Q s is defined by A
sun
jm jrðXþ ; Xsun Þetext ðjrs rTOA jÞ . (6) p sun In Eqs. (1)–(6), sext and sscat are the extinction and the scattering coefficients, respectively, P is the phase function, F sun is the incident solar flux, and A and r are the surface albedo and the normalized bi-directional reflection function, respectively. The solar optical depth tsun ext ðjr rTOA jÞ between the generic point r and the characteristic point at the top of the atmosphere rTOA is computed in a spherical atmosphere. The direction X is characterized by the zenith angle y and the azimuthal angle j, while m stands for the cosine of the zenith angle. Our convention is that m ¼ 1 for upwelling radiation and m ¼ 1 for downwelling Q s ðr s ; Xþ Þ ¼ F sun
radiation. An upward direction is denoted by Xþ, while the notation X is used for a downward direction. The boundary-value problem for the radiative transfer can be cast in an operator form. For this purpose we define the forward radiative transfer operator Z dI s ðrÞ LIðr; XÞ ¼ m ðr; XÞ þ sext ðrÞIðr; XÞ scat Pðr; X; X0 ÞIðr; X0 Þ dO0 dr 4p 4p Z A dðr r s ÞHðmÞm Iðr; X0 ÞHðm0 Þjm0 jrðX; X0 Þ dO0 (7)
p
4p
and the forward source function Q ðr; XÞ ¼ Jss ðr; XÞ þ dðr r s ÞHðmÞmQ s ðr; XÞ, where d is the Dirac delta function and H is the Heaviside step function. The operator representation theorem for the radiative transfer states that the diffuse radiation field Iðr; XÞ, solving the inhomogeneous operator equation
LIðr; XÞ ¼ Q ðr; XÞ,
(8)
ARTICLE IN PRESS 466
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
with the boundary conditions Iðr TOA ; X Þ ¼ Iðr s ; Xþ Þ ¼ 0,
(9)
also solves the radiative transfer equation (1) with the boundary conditions (2) and (3), and vice versa. This assertion follows from the equivalence between the boundary-value problems (1)–(3) and (8)–(9), and the integral representations of the radiative transfer equation for upwelling and downwelling radiation (see Appendix A). It is worth to notice that defining the forward radiative transfer operator as in Eq. (7) we essentially transform a differential equation with inhomogeneous boundary conditions into a differential equation with homogeneous boundary conditions.
2.2. The adjoint radiative transfer operator The adjoint radiative transfer operator Ly is defined through the relation hLI; Iy i ¼ hI; Ly Iy i,
(10)
where the scalar product of the fields I1 and I2 is given by Z rTOA Z hI1 ; I2 i ¼ I1 ðr; XÞI2 ðr; XÞ dO dr. rs
4p
The expression of the adjoint operator Ly, under the assumptions that Iðr TOA ; X Þ ¼ Iðr s ; Xþ Þ ¼ 0, and that Iy ðr TOA ; Xþ Þ ¼ Iy ðr s ; X Þ ¼ 0 is given by Z dIy sscat ðrÞ ðr; XÞ þ sext ðrÞIy ðr; XÞ Pðr; X; X0 ÞIy ðr; X0 Þ dO0 4p dr 4p Z A dðr r s ÞHðmÞjmj Iy ðr; X0 ÞHðm0 Þjm0 jrðX; X0 Þ dO0 .
Ly Iy ðr; XÞ ¼ m
p
(11)
4p
The constructive derivation of Eq. (11) relies on the representation of the forward radiative transfer operator as a linear combination of four operators (cf. Eq. (7)), and on the application of the Lagrange identity (10). The requirement of homogeneous boundary conditions is essential when computing the adjoint of the differential operator of the first order md=dr by using the integration by parts rule.
2.3. The adjoint radiative transfer equation In this section we formulate the adjoint radiative transfer equation in order to express the radiance measured by a satellite instrument as the scalar product of the solution of the adjoint radiative transfer equation and the forward source function. The specification of the adjoint source function and implicitly, the formulation of the adjoint radiative transfer equation, is performed separately for a nadir and a limb viewing geometry [9]. Let us consider a nadir viewing geometry with the line of sight bounded by the surface point S and the point at the top of the atmosphere A. The diffuse radiance at the top of the atmosphere r A ¼ r TOA , and in the measurement direction Xm can be computed by integrating the radiative transfer equation along the line of sight, Z IðrA ; Xm Þ ¼ Iðrs ; Xm Þetext ðjrA rs jÞ þ Jðr; Xm Þetext ðjrA rjÞ ds, jrA rs j
where Iðrs ; Xm Þ is given by Eq. (3) and Z text ðjr1 r2 jÞ ¼ sext ðr0 Þ ds0 jr1 r2 j
is the extinction optical depth between the points r1 and r2 . In a spherical atmosphere the local solar zenith angle varies along the line of sight and as a result, the computation of the radiance field Iðr; XÞ requires the solution of several plane-parallel problems. An efficient and accurate radiance calculation relies on the separation of the single and the multiple scattering contributions, and on the computation of the multiple scattering radiance for a planeparallel atmosphere and a reference solar zenith angle. For multiple scattering calculation, we essentially employ the
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
467
integral approximation Z
f ðr; Xm ; XÞIðr; XÞds ¼
Z
r2
f ðr; Xm ; XÞIðr; XÞ
r1
jr1 r2 j
1
mr
Z
dr
mðrÞ
r2
f ðr; Xmr ; XÞIðr; XÞ dr, r1
where f is an arbitrary scalar function, X is an arbitrary scattering direction and mr ¼ cos yr . The radiance Iðr; XÞ corresponds to a reference radial line, while the overly zealous notation Xmr ¼ ðyr ; jr Þ has been introduced in order to evidence that the polar angles of the measurement direction Xm are specified in a local spherical coordinate system along the reference radial line. In principle, the reference radial line may correspond to any point on the line of sight, but numerical experiments have shown that the selection of the surface point yields the most accurate results. The nadir radiance can then be expressed as Z r etext ðrTOA Þ Iðr s ; X Þjm jrðXmr ; X Þ dO p 2p Z Z rTOA 1 sscat ðrÞ þ Pðr; Xmr ; X0 ÞIðr; X0 Þ dO0 4p mr rs 4p A
IðrA ; Xm Þ ¼
r
r
e½text ðrTOA Þtext ðrÞ dr þ Iss ðrA ; Xm Þ,
(12)
where Iss ðrA ; Xm Þ ¼ Q s ðr s ; Xm Þetext ðjrA rs jÞ þ
Z
J ss ðr; Xm Þetext ðjrA rjÞ ds jrA rs j
is the single scattering radiance and
trext ðrÞ ¼
Z
1
mr
r
sext ðrÞ dr
rs
is the slant extinction optical depth along the reference direction. The above representation shows that the single scattering radiance is calculated for a full spherical geometry, while the multiple scattering radiance is computed for a plane-parallel geometry. Defining the adjoint source function for the nadir radiance by the relation Q y ðr; XÞ ¼
1 sscat ðrÞ A r r r Pðr; Xmr ; XÞe½text ðrTOA Þtext ðrÞ þ dðr r s ÞHðmÞjmjrðXmr ; XÞetext ðrTOA Þ , 4p p
mr
the integral representation Im ¼ hQ y ; Ii þ Iss ,
(13)
with Im ¼ IðrA ; Xm Þ, holds true. It should be pointed out that the equivalence between the representations (13) and (12) follows directly from the definition of the scalar product hQ y ; Ii ¼
Z
Z
rTOA
Q y ðr; XÞIðr; XÞ dO dr,
4p
rs
and the sifting property of the Dirac delta function, which yields A
p
r
etext ðrTOA Þ
Z
rTOA
dðr rs Þ dr rs
¼
A
p
r
etext ðrTOA Þ
Z
Iðr; XÞHðmÞjmjrðXmr ; XÞ dO 4p
Z
Iðr s ; X Þjm jrðXmr ; X Þ dO . 2p
Let us now consider a limb viewing geometry with the line of sight bounded by the near side point A and the far side point A0 , and let the radial coordinate of the tangent point of the line of sight T be r T . The radiance measured by a limb instrument is given by IðrA ; Xm Þ ¼
Z
sscat ðrÞ 4p jrA rA0 j
Z 4p
Pðr; Xm ; X0 ÞIðr; X0 Þ dO0 etext ðjrA rjÞ ds þ Iss ðrA ; Xm Þ,
where Iss ðrA ; Xm Þ ¼
Z
J ss ðr; Xm Þetext ðjrA rjÞ ds. jrA rA0 j
(14)
ARTICLE IN PRESS 468
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
m
A M
r T
r+
M’
rTOA rT
r-
A’
rTOA
Fig. 1. Limb scattering geometry.
The line integral in Eq. (14) can be transformed into a simple integral along the radial coordinate by using the basic result Z Z Z f ðrÞ ds ¼ f ðrþ Þ ds þ f ðr Þ ds jrA rA0 j
jrA rT j rA
Z ¼
rT
jr r j
T A ds ½f ðr Þ þ f ðr Þ ðrÞ dr, dr þ
where rþ and r , with jrþ j ¼ jr j ¼ r, are ‘‘symmetric’’ integration points situated on the near side and on the far side of the line of sight, respectively (Fig. 1), and ds 1 r ðrÞ ¼ ¼ Hðr r T Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . dr mðrÞ 2 r r2 T
As for a nadir viewing geometry, the multiple scattering radiance is computed in a plane-parallel atmosphere and for a reference solar zenith angle. The variation with the solar zenith angle is accounted for in an approximate manner by scaling the reference radiance, that is, Pðr; Xm ; XÞIðr; XÞ cðrÞPðr; Xmr ; XÞIðr; XÞ.
(15)
The scaling factor is computed as sun
cðrÞ ¼
etext ðjrrTOA jÞ , sun etext ðjrr rTOA jÞ
where rr is a vector along the radial line for which it holds true that jrr j ¼ jrj ¼ r. The exponential term expftsun ext ðjr rTOA jÞg reproduces the dependency of the source term on the solar zenith angle and is assumed to act like a multiplicative factor for the radiance fields. Noticing that in practice, the reference calculation corresponds to the tangent point, we obtain Z Z rA ds sscat ðrÞ IðrA ; Xm Þ ¼ Pðr; Xmr ; XÞIðr; XÞ dO0 TðrÞ ðrÞ dr þ Iss ðrA ; Xm Þ, (16) 4p dr rT 4p with TðrÞ ¼ cðrþ Þetext ðjrA r
þ
jÞ
þ cðr Þetext ðjrA r jÞ .
In view of Eq. (16), the adjoint source function for the limb radiance is defined by ds sscat ðrÞ Pðr; Xmr ; XÞTðrÞ ðrÞ, Q y ðr; XÞ ¼ Hðr r T Þ 4p dr in which case, IðrA ; Xm Þ possesses the integral representation (13). We are now in the position to formulate the main result of the adjoint radiative transfer approach: If Iðr; XÞ solves the forward problem
LIðr; XÞ ¼ Q ðr; XÞ,
(17)
Iðr TOA ; X Þ ¼ Iðr s ; Xþ Þ ¼ 0,
(18)
y
and I ðr; XÞ solves the adjoint problem
Ly Iy ðr; XÞ ¼ Q y ðr; XÞ,
(19)
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Iy ðr TOA ; Xþ Þ ¼ Iy ðr s ; X Þ ¼ 0,
469
(20)
then the radiance measured by a satellite instrument can be expressed as ð13Þ
Im ¼ hQ y ; Ii þ Iss ð19Þ
¼ hLy Iy ; Ii þ Iss
ð10Þ;ð18Þ;ð20Þ
¼
hIy ; LIi þ Iss
ð17Þ
¼ hIy ; Q i þ Iss .
(21)
The scalar product representation (21) reveals that the information about the measurement geometry is encapsulated in the adjoint source function Q y (and implicitly in the adjoint radiation field Iy ), while the information about the source conditions are contained in the forward source function Q . 3. Solution of the adjoint radiative transfer equation The forward and the adjoint radiative transfer equations (17) and (19) are related to each other. Replacing X by X in the expression of the adjoint operator Ly, we rewrite the adjoint radiative transfer equation as
Ly ðXÞIy ðr; XÞ ¼ Q y ðr; XÞ. y Defining the conjugate adjoint field ^I by the relation
^Iy ðr; XÞ ¼ Iy ðr; XÞ
(22) 0
0
and using the symmetry property of the phase function Pðr; X; X Þ ¼ Pðr; X ; XÞ, we find the identity [6] y
Ly ðXÞIy ðr; XÞ ¼ L^I ðr; XÞ.
(23)
The solution of the adjoint radiative transfer equation can be found by using the same solution method as for the forward problem with an appropriate source function and boundary conditions. Performing an azimuthal expansion of the radiance field we obtain the boundary-value problems:
m
y Z y y d^Im sscat ðrÞ 1 ðr; mÞ ¼ sext ðrÞ^Im ðr; mÞ þ pm ðr; m; m0 Þ^Im ðr; m0 Þ dm0 2 dr 1 sscat ðrÞ r r þ ð2 dm0 ÞF r pm ðr; m; mr Þe½text ðrTOA Þtext ðrÞ , 4p
^Iy ðr TOA ; m Þ ¼ 0, m ^Iy ðr s ; mþ Þ ¼ 2A m
Z
0
1
^Iy ðr s ; m Þjm jr ðmþ ; m Þ dm þ ð2 dm0 ÞF r A m r ðmþ ; m Þetrext ðrTOA Þ , r m r m m
p
(24)
for a nadir viewing geometry, and
m
y Z y y dI^m sscat ðrÞ 1 ðr; mÞ ¼ sext ðrÞ^Im ðr; mÞ þ pm ðr; m; m0 Þ^Im ðr; m0 Þ dm0 2 dr 1 ds sscat ðrÞ þ ð2 dm0 ÞHðr r T Þ pm ðr; m; mr ÞTðrÞ ðrÞ, 4p dr
^Iy ðr TOA ; m Þ ¼ 0, m ^Iy ðr s ; mþ Þ ¼ 2A m
Z
0
1
^Iy ðr s ; m Þjm jr ðmþ ; m Þ dm m m
(25)
for a limb viewing geometry. In Eqs. (24) and (25), pm and rm are the azimuthal expansion coefficients of the phase function and the bi-directional reflection function, respectively, and F r ¼ 1=mr . It should be pointed out that the definition of the y conjugate adjoint field (22) implies that Iym ðr; mÞ ¼ ^Im ðr; mÞ. The adjoint radiative transfer problems (24) and (25) are solved by using the discrete ordinate method with matrix exponential [13]. The peculiarities of this solution method for the present application will be sketched below. The atmosphere is discretized in N levels: r 1 4r 2 4 4r N with the convention r 1 ¼ r TOA and r N ¼ r s . A layer j, bounded above by the level r j and below by the level r jþ1 , has the geometrical thickness Dr j ¼ r j r jþ1 . The extinction and the scattering coefficients as well as the phase function coefficients are assumed to be constant within each layer; their average values are s¯ ext;j , s¯ scat;j and p¯ m;j , respectively. The angular variation of the phase function and radiance is discretized in M
ARTICLE IN PRESS 470
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
discrete ordinate directions per hemisphere, and the set of Gauss–Legendre quadrature points and weights in the interval ð0; 1Þ is denoted by fmk ; wk gk¼1;M . The matrix exponential formalism operates with the concept of layer equation
A1m;j im;j þ A2m;j im;jþ1 ¼ bm;j ,
(26)
which relates the level values of the adjoint radiance field, 2 1 3 im;j im;j ¼ 4 2 5, im;j 1;2
y
with ½im;j k ¼ ^Im ðr j ; mk Þ and k ¼ 1; . . . ; M. The layer matrices can be expressed in compact form as
A1m;j ¼ D1m;j V1 m;j , A2m;j ¼ D2m;j V1 m;j , where the columns of the matrix Vm;j are the eigenvectors of the layer matrix 2 11 3 Am;j A12 m;j 4 5, Am;j ¼ A12 A11 m;j m;j with entries 1 p¯ ðm ; m Þ 2s¯ ext;j dkl , ½w s¯ 2mk l scat;j m;j k l 1 ½A12 p¯ ðm ; ml Þ, w s¯ m;j kl ¼ 2mk l scat;j m;j k
½A11 m;j kl ¼
and the diagonal scaling matrices D1m;j and D2m;j are given by D1m;j ¼ diag½a0 ðlk Dr j Þ; 1, D2m;j ¼ diag½1; a0 ðlk Dr j Þ, with lk , k ¼ 1; . . . ; M, being the positive eigenvalues of the matrix Am;j and a0 ðxÞ ¼ ex . The expression of the layer vector depends on the viewing geometry. For nadir measurements, bm;j is computed as
bm;j ¼ Lm;j V1 m;j ðDr j bm;j Þ, where the components of the source vector 2 1 3 bm;j bm;j ¼ 4 2 5 bm;j are given by 1;2
½bm;j k ¼
1
mk
ð2 dm0 Þ
Fr p¯ ðmk ; mr Þ. s¯ 4p scat;j m;j
The diagonal scaling matrix Lm;j is expressed as
Lm;j ¼ diag½b1 ðlk Dr j Þ; b2 ðlk Dr j Þ,
(27)
with the interpolation functions b1 and b2 being defined by r
b1 ðxÞ ¼ b2 ðxÞ ¼
r
eðtext;j þxÞ etext;jþ1 , trext;jþ1 trext;j x ðtrext;jþ1 þxÞ
e
trext;j
e
trext;j trext;jþ1 x
.
(28)
The layer equation has been derived by assuming that the extinction optical depth varies linearly within the layers; the level values of the slant extinction optical depth being given by
trext;j ¼
j1 1X
mr u¼1
s¯ ext;u Dru ;
j ¼ 1; . . . ; N.
For a limb viewing geometry, we consider a division of the line of sight in an even number N point LOS of points; the number of point point ¯ limb paths along the line of sight is then Npath LOS ¼ N LOS 1. The index j ¼ N LOS =2 þ 1 corresponds to the tangent point of radial coordinate r T , and for simplicity, we assume that ¯j is also the level index for which it holds true that r ¯j ¼ r T . The
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
471
¯ number of effective layers traversed by the line of sight is Neff lay ¼ j 1, and it is apparent that bm;j vanishes for all layers j4Neff lay . To derive the expression of the layer vector we introduce the following quantities:
(1) the scaling factors sun
cj ¼
etext ðjrj rTOA jÞ ; sun etext ðjrr;j rTOA jÞ
j ¼ 1; . . . ; Npoint LOS ,
where jrj j ¼ jrr;j j ¼ r levðjÞ , and the integer function levðjÞ maps the point index j, into the level index ( j; 1pjp¯j; levðjÞ ¼ ¯ ¯ 2j j; j þ 1pjp2¯j 1; (2) the slant optical depths
trext;1 ¼ 0;
trext;j ¼
j1 X
s¯ ext;layðuÞ Du ; j ¼ 2; . . . ; Npoint LOS ,
u¼1
where Dj are the limb path lengths, and the integer function layðuÞ map the path index u, into the layer index ( u; 1pup¯j 1; layðuÞ ¼ 2¯j u 1; ¯jpup2¯j 2; (3) the average Jacobian values
¯ ds dr
¼ j
Dj ; DrlayðjÞ
j ¼ 1; . . . ; Npath LOS .
The layer vector entering in Eq. (26) can then be expressed as þ 1 1 bm;j ¼ Lþ m;j V m;j ðDr j bm;j Þ þ Lm;j Vm;j ðDr j bm;j Þ, þ
where the components of the source vectors bm;j and bm;j , are given by þ1;2
½bm;j k ¼
1
mk
ð2 dm0 Þ
¯ 1 ds s¯ scat;j p¯ m;j ðmk ; mr Þ¯cj 4p dr j
ð2 dm0 Þ
¯ 1 ds s¯ scat;j p¯ m;j ðmk ; mr Þ¯c2¯jj1 , 4p dr 2¯jj1
and 1;2
½bm;j k ¼
1
mk
respectively. The diagonal matrix Lþ m;j is as in Eqs. (27) and (28), while the interpolation functions b1 and b2 , which enter in r r the expression of the diagonal matrix L m;j , are expressed in terms of the slant optical depths text;2¯jj and text;2¯jj1 in place r r of text;j and text;jþ1 , respectively. The layer equation (26) together with the boundary conditions at the top and the bottom of the atmosphere are assembled into the global matrix of the entire atmosphere, and the solution of the resulting system of equations yields the level values of the adjoint radiance field. An essential aspect of the adjoint approach is the calculation of the scalar product in Eq. (21). The integration with respect to the radial coordinate can be analytically performed (see Appendix B) and this result substantially increases the efficiency of the method. 4. Numerical simulations The purpose of our numerical simulations is to validate the scalar product representation (21), since this equation is the basic relation for computing the weighting functions by the adjoint approach. The radiance measured by a satellite instrument is computed by solving the forward and the adjoint problems in the framework of the discrete ordinate method with matrix exponential. The calculations are performed in the spectral domain between 320 and 350 nm, and for a vertically inhomogeneous atmosphere taking into account Rayleigh and aerosol scattering, as well as the absorption by O3, NO2 and aerosol. The O3 and NO2 profiles have been taken from the model of McLinden et al. [14], while the LOWTRAN/MODTRAN database has been used for aerosol parametrization [15]. The atmosphere is discretized with a step of 1 km between 0 and 25 km, a step of 2.5 km between 25 and 50 km and a step of 5 km between 50 and 100 km. The surface albedo is set to 0.3, and the solar flux is 1.
ARTICLE IN PRESS 472
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Radiance [relative units]
2e−01
1 2 3 4 5 6 7 8
1e−01
0e+00
1 2 3 4 5 6 7 8
320
340
330
350
Relative Difference
2e−03
8 1e−03 8
7 4 3 5 2 1 6
12 3 4 5 6
7
0e+00 320
330
340
350
Wavelength [nm] Fig. 2. (a) Nadir radiances for different values of the solar zenith angle. (b) Relative errors between the solutions computed by using the forward and the adjoint method. The curves correspond to the following values of the solar zenith angle: 1–10 , 2–20 , 3–30 , 4–40 , 5–50 , 6–60 , 7–70 and 8–80 . The local viewing zenith angle is 10 and all angles refer to the top of the atmosphere.
The results shown in Fig. 2 correspond to a nadir viewing geometry. In Fig. 2a we plot the radiances computed by the forward method for different values of the solar zenith angle. The relative errors between the forward and the adjoint methods are smaller than 0.13% (Fig. 2b), while the computational costs are similar. In Fig. 3 we illustrate the simulation results for a limb viewing geometry. For the forward method, the radiative transfer equation is solved for five values of the solar zenith angle (five right-hand side vectors of the matrix equation), while, in view of the approximation (15), a single calculation is performed for the adjoint method. Due to the approximation employed, the relative errors are of about 0.8%. The computer time needed to simulate the radiances at 121 wavelengths and eight discrete ordinate directions per hemisphere was 24 s for the forward method and 14 s for the adjoint method (on an Intel Xeon processor with 2.8 GHz). The main problem of the adjoint method is that for large optical depths along the line of sight when the tangent point is close to the terminator, the approximation (15) deteriorates and the relative errors may reach 3% (Fig. 4). However, in atmospheric remote sensing, due to the general nature of the Gauss–Newton iterative method commonly used to solve the nonlinear inverse problem, accuracies of the weighting functions up to a few percent do not deteriorate the convergence rate of the solution, and the application of the adjoint approach is not critical.
5. Conclusions The adjoint radiative transfer is formulated for a pseudo-spherical atmosphere and general viewing geometries. The single scattering radiance is computed in a spherical atmosphere, while for the multiple scattering radiance we formulate an adjoint radiative transfer equation in a plane-parallel atmosphere starting from the integral representation of the radiance measured by a satellite instrument. For limb scattering geometries, we solve the one-dimensional radiative transfer problem for a reference solar zenith angle, and scale the calculated radiances by a correction factor depending on the solar optical depth. The solution of the adjoint radiative transfer equation is obtained by using the discrete ordinate method with matrix exponential. The accuracy of the adjoint model is satisfactory for a nadir viewing geometry and still acceptable for a limb viewing geometry. The adjoint approach together with a linearized forward model are implemented in a software package that we developed at the German Aerospace Centre for solving inverse problems in atmospheric remote sensing. A comparison of these two linearization methods with respect to their efficiency and accuracy will be the subject of a future paper.
ARTICLE IN PRESS
Radiance [relative units]
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
473
2e−01
96 5 78
1e−01
4 3
9 8 7 6 5 4 3
0e+00 12 320
12
330
340
350
340
350
Relative Difference
8e−03 9 8 7 4 6e−03 56
9 8 6 45 7
4e−03 3
3 12
2e−03 2 1
0e−00 320
330 Wavelength [nm]
Fig. 3. (a) Limb radiances (b) Relative errors between the solutions computed by using the forward and the adjoint method. The curves correspond to the following values of the tangent height: 1–50, 2–40, 3–30, 4–20; 5–16; 6–13; 7–10; 8–6; 9–3 km. The solar zenith angle is 45 and the relative azimuthal angle is 45 .
Radiance [relative units]
2e−01 98765
1e−01 4
0e+00
3
5 7 6 9 8 4 3 2 1
320
21
330
340
350
340
350
Relative Difference
3e−02 9 8
2e−02
9 8
7 6 5
7
1e−02
4
6
3 1 2
5 12 3 4
0e+00 320
330 Wavelength [nm]
Fig. 4. Same as in Fig. 3 but for a solar zenith angle of 80 and a relative azimuthal angle of 0 .
ARTICLE IN PRESS 474
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Acknowledgement The authors are grateful to Eugene Ustinov for illuminating discussions regarding the modeling of limb radiances in the framework of adjoint radiative transfer. Appendix A In this appendix we prove that a diffuse radiation field solving the inhomogeneous operator equation (8) with the boundary conditions (9) also solves the radiative transfer equation (1) with the boundary conditions (2) and (3). The proof relies on the equivalence between the boundary-value problem (1)–(3) and the integral representations of the radiative transfer equation for upwelling and downwelling radiation: If the diffuse radiation field Iðr; XÞ satisfies the boundary-value problem (1)–(3), then Iðr; XÞ solves the integral equation Z rTOA 1 0 Iðr; X Þ ¼ Jðr 0 ; X Þe½text ðr Þtext ðrÞ dr 0 (29) jm j r for X ¼ X, with
text ðrÞ ¼
Z
1 jm
j
r
sext ðrÞ dr rs
and þ
Iðr; Xþ Þ ¼ Iðr s ; Xþ Þetext ðrÞ þ
1
m
þ
Z
r
þ
þ
Jðr 0 ; Xþ Þe½text ðrÞtext ðr Þ dr 0 0
(30)
rs
for X ¼ Xþ, with
tþext ðrÞ ¼
1
m
þ
Z
r
sext ðrÞ dr. rs
The converse results is also true. Let us now assume that the diffuse radiation field Iðr; XÞ satisfies the boundary-value problem (8)–(9). For downwelling radiances, there hold X ¼ X and m ¼ m o0. As a result we deduce that HðmÞ ¼ 0, and the differential equation (8) can be written explicitly as
m
dI ðr; X Þ ¼ sext ðrÞIðr; X Þ þ Jðr; X Þ. dr
(31)
Integrating the inhomogeneous differential equation (31) with the boundary condition Iðr TOA ; X Þ ¼ 0; yields the integral form of the radiative transfer equation for the downwelling radiation (29). For upwelling radiances, the conditions X ¼ Xþ and m ¼ mþ 40 imply that HðmÞ ¼ 1, and the differential equation (8) read as Z dI A mþ ðr; Xþ Þ ¼ sext ðrÞIðr; Xþ Þ þ Jðr; Xþ Þ þ dðr rs Þmþ Iðr; X Þjm jrðXþ ; X Þ dO dr p 2p þ dðr r s Þmþ Q s ðr; Xþ Þ.
(32) þ
Integrating the inhomogeneous differential equation (32) with the boundary condition Iðr s ; X Þ ¼ 0, and using the results Z r Z þ þ 1 A 0 dðr 0 rs Þmþ Iðr 0 ; X Þjm jrðXþ ; X Þ dO e½text ðrÞtext ðr Þ dr 0 þ m rs p 2p Z þ A ¼ Iðr s ; X Þjm jrðXþ ; X Þ dO etext ðrÞ
p
2p
and 1
mþ
Z
r rs
þ
þ
þ
½dðr 0 r s Þmþ Q s ðr 0 ; Xþ Þe½text ðrÞtext ðr Þ dr 0 ¼ Q s ðr s ; Xþ Þetext ðrÞ , 0
we find the integral form of the radiative transfer equation for the upwelling radiation (30). Appendix B In this appendix we sketch the calculation of the scalar product in Eq. (21) for a nadir viewing geometry; the treatment for a limb viewing geometry is similarly. Considering the multiple scattering radiance Ims ¼ hIy ; Q i, and assuming the
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
475
azimuthal expansions 2M1 X Ims ðXr Þ ¼ Ims;m ðmr Þ cos mjr , m¼0
Iy ðr; XÞ ¼
2M1 X
Iym ðr; mÞ cos½mðj jr Þ,
m¼0
Q ðr; XÞ ¼
2M1 X
Q m ðr; mÞ cos mj,
m¼0
we find that Ims;m ¼ ð1 þ dm0 Þp
Z
r TOA rs
Z
1
Q m ðr; mÞIym ðr; mÞ dm dr.
1
Using the expression of the mth Fourier harmonics of the source function, sscat ðrÞ sun Q m ðr; mÞ ¼ ð2 dm0 ÞF sun pm ðr; m; msun Þetext ðjrrTOA jÞ 4p A sun þ ð2 dm0 ÞF sun dðr r s ÞHðmÞmjmsun jrm ðm; msun Þetext ðjrs rTOA jÞ ,
p
we obtain 1 sun 1 NX 2 F sun sTm;j Im;j þ 2AF sun jmsun jetext;N rTm im;N , 2 j¼1
Ims;m ¼
where the components of the scattering vector sm;j are given by ¼ wl s¯ scat;j p¯ m;j ð ml ; msun Þ, ½s1;2 m;j l the reflection vector rm has the entries ½rm l ¼ wl ml rm ðml ; msun Þ, and accordingly to the average secant approximation, Im;j is the integral Z Dr j ½ðr=Dr j Þtsun þð1r=Drj Þtsun ext;j ext;jþ1 i Im;j ¼ e (33) m;j ðrÞ dr. 0
To compute Im;j we use the analytical representation of the internal radiance derived in [13], that is, 2 3 " # 1 diag½elk ðDrj rÞ xm;j diag½B1 ðlk ðDr j rÞÞ 0 5 þ Vm;j V1 Dr b , im ðrÞ ¼ Vm;j 4 2 0 diag½B2 ðlk rÞ m;j j m;j diag½elk r xm;jþ1 1
2
1
2
1 T with xm;j ¼ ½xm;j ; xm;j T ¼ V1 m;j im;j , xm;jþ1 ¼ ½xm;jþ1 ; xm;jþ1 ¼ V m;j im;jþ1 , and ðtrext;j þxÞ
B1 ðxÞ ¼
e
½ðr=Dr j Þtrext;j þð1r=Dr j Þtrext;jþ1
e
r
B2 ðxÞ ¼
,
trext;jþ1 trext;j lk Drj r
r
eðtext;jþ1 þxÞ e½ðr=Drj Þtext;j þð1r=Drj Þtext;jþ1 . trext;j trext;jþ1 lk Drj
The result of the integration is 2 3 " 1 ½diagða1 ðlk Dr j ÞÞxm;j ½diagðB¯ 1 ðlk Dr j ÞÞ 5 þ Vm;j Im;j ¼ Vm;j 4 2 0 ½diagða2 ðlk Dr j ÞÞxm;jþ1
0 ½diagðB¯ 2 ðlk Dr j ÞÞ
# V1 m;j Dr j bm;j ,
where sun a1 ðxÞ ¼ Dr j I0 ðtsun ext;j ; text;jþ1 þ xÞ, sun a2 ðxÞ ¼ Dr j I0 ðtsun ext;j þ x; text;jþ1 Þ,
B¯ 1 ðxÞ ¼ B¯ 2 ðxÞ ¼
Drj
trext;jþ1 trext;j lk Dr j Drj
trext;j trext;jþ1 lk Dr j
trext;j
½e
r
sun sun r sun r ½etext;jþ1 I0 ðtsun ext;j þ x; text;jþ1 Þ I0 ðtext;j þ text;j ; text;jþ1 þ text;jþ1 Þ,
and I0 is the integral Z 1 I0 ða; bÞ ¼ e½axþbð1xÞ dx ¼ 0
sun sun r sun r I0 ðtsun ext;j ; text;jþ1 þ xÞ I0 ðtext;j þ text;j ; text;jþ1 þ text;jþ1 Þ,
1 ðea eb Þ. ab
References [1] Marchuk GI. Equation for the value of information from weather satellites and formulation of inverse problems. Cosmic Res 1964;2: 394–409.
ARTICLE IN PRESS 476
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
[2] Box MA, Gerstl SAW, Simmer C. Application of the adjoint formulation to the calculation of atmospheric radiative effects. Beitr Phys Atmos 1988;61:303–11. [3] Ustinov EA. Inverse problem of photometric observation of solar radiation reflected by an optically dense planetary atmosphere: mathematical methods and weighting functions of linearized inverse problem. Cosmic Res 1991;29:519–32. [4] Ustinov EA. Adjoint sensitivity analysis of radiative transfer equation: temperature and gas mixing ratio weighting functions for remote sensing of scattering atmospheres in thermal IR. JQSRT 2001;68:195–211. [5] Ustinov EA. Atmospheric weighting functions and surface partial derivatives for remote sensing of scattering planetary atmospheres in thermal spectral region: general adjoint approach. JQSRT 2005;92:351–71. [6] Rozanov VV, Rozanov AV. Generalized form of the direct and adjoint radiative transfer equations. JQSRT 2007;104:155–70. [7] Landgraf J, Hasekamp OP, Box MA, Trautmann T. A linearized radiative transfer model for ozone profile retrieval using the analytical forward-adjoint perturbation theory approach. J Geophys Res 2001;106(D21):27291–305. [8] Walter H, Landgraf J, Hasekamp OP. Linearization of a pseudo-spherical vector radiative transfer model. JQSRT 2004;85:251–83. [9] Ustinov EA. Adjoint approach to evaluation of weighting functions for remote sensing of scattering planetary atmospheres in thermal spectral region with limb-viewing geometry. EGU General Assembly 2008, Vienna, Geophysical Research Abstracts, vol. 10, EGU2008-A-10664, 2008, SRef-ID: 16077962/gra/EGU2008-A-10664. [10] Rozanov VV, Kurosu T, Burrows JP. Retrieval of atmospheric constituents in the UV-visible: A new quasi-analytical approach for the calculation of the weighting functions. JQSRT 1998;60:277–99. [11] Rozanov VV, Kokhanovsky AA. The average number of photon scattering events in vertically inhomogeneous atmosphere. JQSRT 2005;96:11–33. [12] Qin Y, Jupp DLB, Box MA. Extension of the discrete-ordinate algorithm and efficient radiative transfer calculation. JQSRT 2002;74:767–81. [13] Doicu A, Trautmann T. Discrete ordinate method with matrix exponential for a pseudo-spherical atmosphere. Scalar case. JQSRT 2009;110:146–58. [14] McLinden CA, McConnell JC, Griffioen E, McElroy CT. A vector radiative-transfer model for the Odin/OSIRIS project. Can J Phys 2002;80:375–93. [15] Kneizys F, Abreu L, Anderson G, Chetwynd J, Shettle E, Berk A, Bernstein L, Robertson D, Acharya P, Rothman L, Selby J, Gallery W, Clough S. The MODTRAN 2/3 report and LOWTRAN 7 model. Technical Report, contract F19628-91-C-0132 with Ontar Corporation, Phillips Laboratory, Hanscom AFB; 1996.
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
Adjoint problem of radiative transfer for a pseudo-spherical atmosphere and general viewing geometries A. Doicu , T. Trautmann Remote Sensing Technology Institute, German Aerospace Centre, Oberpfaffenhofen, Wessling, Germany
a r t i c l e i n f o
abstract
Article history: Received 14 October 2008 Received in revised form 19 January 2009 Accepted 19 January 2009
We formulate the adjoint radiative transfer for a pseudo-spherical atmosphere and various retrieval scenarios. The single scattering radiance is computed in a spherical atmosphere by using the source integration technique, while for the multiple scattering radiance we formulate an one-dimensional adjoint radiative transfer equation in a plane-parallel atmosphere. The adjoint solution of the radiative transfer equation is obtained by employing the discrete ordinate method with matrix exponential. We provide an abbreviated derivation of our formalism as well as a discussion of the numerical implementation of the theory. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Radiative transfer Pseudo-spherical atmosphere Adjoint operator
1. Introduction The adjoint approach allows the linearization of radiative transfer problem with respect to the optical properties of the atmosphere and the surface. The application of this approach to weighting functions calculation requires the formulation and the solution of the adjoint radiative transfer equation. For a plane-parallel atmosphere, formulations of the adjoint radiative transfer have been given by Marchuk [1], Box et al. [2], Ustinov [3–5], Rozanov and Rozanov [6] and Landgraf et al. [7]. For a pseudo-spherical atmosphere, Walter et al. [8] applied the adjoint approach to a nadir viewing geometry, in which case, the singly scattered radiation is calculated for a spherical geometry, while the multiply scattered radiation is computed for a plane-parallel geometry. Recently, Ustinov [9] extended the general adjoint approach for nadir geometries to limb geometries, but determined the total scattered radiation in a plane-parallel manner. Several solution methods have been employed for solving the adjoint radiative transfer problem. In this context we mention the finite difference method [10], the discrete ordinate method [11,12], and the Gauss–Seidel iteration technique [7]. The goal of our paper is to formulate the adjoint radiative transfer for a pseudo-spherical atmosphere, and for both nadir and limb viewing geometries. The adopted solution method for solving the adjoint radiative transfer problem is the discrete ordinate method with matrix exponential [13]. This method is unconditionally stable for arbitrarily large optical depths and can be simply extended to the present application.
Corresponding author.
E-mail address: [email protected] (A. Doicu). 0022-4073/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2009.01.027
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
465
2. The adjoint radiative transfer approach Depending on the choice of the forward operator and the treatment of boundary conditions, equivalent formulations of the adjoint radiative transfer can be derived [4,6]. In the present analysis we adopt the technique described in [4,7], which assumes the following steps: (1) representation of the boundary-value problem for the radiative transfer in an operator form; (2) derivation of the adjoint radiative transfer operator from the Lagrange identity; (3) derivation of the adjoint radiative transfer problem by using the integral representation of the radiance measured by a satellite instrument. The main results of this section will be written with cursive letters, and for reasons of conciseness, their proofs will be only sketched. 2.1. The forward radiative transfer equation In a pseudo-spherical atmosphere, the boundary-value problem for the diffuse radiance consists in the inhomogeneous differential equation
m
dI ðr; XÞ ¼ sext ðrÞIðr; XÞ þ Jðr; XÞ dr
(1)
the top-of-atmosphere boundary condition (r ¼ r TOA ), Iðr TOA ; X Þ ¼ 0
(2)
and the surface boundary condition (r ¼ r s ), Z A Iðr s ; Xþ Þ ¼ Iðr s ; X Þjm jrðXþ ; X Þ dO þ Q s ðr s ; Xþ Þ.
p
(3)
2p
The source function Jðr; XÞ ¼ J ss ðr; XÞ þ J ms ðr; XÞ includes the contribution of the single and the multiple scattering terms Jss ðr; XÞ ¼ F sun
sscat ðrÞ sun Pðr; X; Xsun Þetext ðjrrTOA jÞ 4p
(4)
and Jms ðr; XÞ ¼
sscat ðrÞ 4p
Z
Pðr; X; X0 ÞIðr; X0 Þ dO0 ,
(5)
4p
respectively, while the surface function Q s is defined by A
sun
jm jrðXþ ; Xsun Þetext ðjrs rTOA jÞ . (6) p sun In Eqs. (1)–(6), sext and sscat are the extinction and the scattering coefficients, respectively, P is the phase function, F sun is the incident solar flux, and A and r are the surface albedo and the normalized bi-directional reflection function, respectively. The solar optical depth tsun ext ðjr rTOA jÞ between the generic point r and the characteristic point at the top of the atmosphere rTOA is computed in a spherical atmosphere. The direction X is characterized by the zenith angle y and the azimuthal angle j, while m stands for the cosine of the zenith angle. Our convention is that m ¼ 1 for upwelling radiation and m ¼ 1 for downwelling Q s ðr s ; Xþ Þ ¼ F sun
radiation. An upward direction is denoted by Xþ, while the notation X is used for a downward direction. The boundary-value problem for the radiative transfer can be cast in an operator form. For this purpose we define the forward radiative transfer operator Z dI s ðrÞ LIðr; XÞ ¼ m ðr; XÞ þ sext ðrÞIðr; XÞ scat Pðr; X; X0 ÞIðr; X0 Þ dO0 dr 4p 4p Z A dðr r s ÞHðmÞm Iðr; X0 ÞHðm0 Þjm0 jrðX; X0 Þ dO0 (7)
p
4p
and the forward source function Q ðr; XÞ ¼ Jss ðr; XÞ þ dðr r s ÞHðmÞmQ s ðr; XÞ, where d is the Dirac delta function and H is the Heaviside step function. The operator representation theorem for the radiative transfer states that the diffuse radiation field Iðr; XÞ, solving the inhomogeneous operator equation
LIðr; XÞ ¼ Q ðr; XÞ,
(8)
ARTICLE IN PRESS 466
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
with the boundary conditions Iðr TOA ; X Þ ¼ Iðr s ; Xþ Þ ¼ 0,
(9)
also solves the radiative transfer equation (1) with the boundary conditions (2) and (3), and vice versa. This assertion follows from the equivalence between the boundary-value problems (1)–(3) and (8)–(9), and the integral representations of the radiative transfer equation for upwelling and downwelling radiation (see Appendix A). It is worth to notice that defining the forward radiative transfer operator as in Eq. (7) we essentially transform a differential equation with inhomogeneous boundary conditions into a differential equation with homogeneous boundary conditions.
2.2. The adjoint radiative transfer operator The adjoint radiative transfer operator Ly is defined through the relation hLI; Iy i ¼ hI; Ly Iy i,
(10)
where the scalar product of the fields I1 and I2 is given by Z rTOA Z hI1 ; I2 i ¼ I1 ðr; XÞI2 ðr; XÞ dO dr. rs
4p
The expression of the adjoint operator Ly, under the assumptions that Iðr TOA ; X Þ ¼ Iðr s ; Xþ Þ ¼ 0, and that Iy ðr TOA ; Xþ Þ ¼ Iy ðr s ; X Þ ¼ 0 is given by Z dIy sscat ðrÞ ðr; XÞ þ sext ðrÞIy ðr; XÞ Pðr; X; X0 ÞIy ðr; X0 Þ dO0 4p dr 4p Z A dðr r s ÞHðmÞjmj Iy ðr; X0 ÞHðm0 Þjm0 jrðX; X0 Þ dO0 .
Ly Iy ðr; XÞ ¼ m
p
(11)
4p
The constructive derivation of Eq. (11) relies on the representation of the forward radiative transfer operator as a linear combination of four operators (cf. Eq. (7)), and on the application of the Lagrange identity (10). The requirement of homogeneous boundary conditions is essential when computing the adjoint of the differential operator of the first order md=dr by using the integration by parts rule.
2.3. The adjoint radiative transfer equation In this section we formulate the adjoint radiative transfer equation in order to express the radiance measured by a satellite instrument as the scalar product of the solution of the adjoint radiative transfer equation and the forward source function. The specification of the adjoint source function and implicitly, the formulation of the adjoint radiative transfer equation, is performed separately for a nadir and a limb viewing geometry [9]. Let us consider a nadir viewing geometry with the line of sight bounded by the surface point S and the point at the top of the atmosphere A. The diffuse radiance at the top of the atmosphere r A ¼ r TOA , and in the measurement direction Xm can be computed by integrating the radiative transfer equation along the line of sight, Z IðrA ; Xm Þ ¼ Iðrs ; Xm Þetext ðjrA rs jÞ þ Jðr; Xm Þetext ðjrA rjÞ ds, jrA rs j
where Iðrs ; Xm Þ is given by Eq. (3) and Z text ðjr1 r2 jÞ ¼ sext ðr0 Þ ds0 jr1 r2 j
is the extinction optical depth between the points r1 and r2 . In a spherical atmosphere the local solar zenith angle varies along the line of sight and as a result, the computation of the radiance field Iðr; XÞ requires the solution of several plane-parallel problems. An efficient and accurate radiance calculation relies on the separation of the single and the multiple scattering contributions, and on the computation of the multiple scattering radiance for a planeparallel atmosphere and a reference solar zenith angle. For multiple scattering calculation, we essentially employ the
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
467
integral approximation Z
f ðr; Xm ; XÞIðr; XÞds ¼
Z
r2
f ðr; Xm ; XÞIðr; XÞ
r1
jr1 r2 j
1
mr
Z
dr
mðrÞ
r2
f ðr; Xmr ; XÞIðr; XÞ dr, r1
where f is an arbitrary scalar function, X is an arbitrary scattering direction and mr ¼ cos yr . The radiance Iðr; XÞ corresponds to a reference radial line, while the overly zealous notation Xmr ¼ ðyr ; jr Þ has been introduced in order to evidence that the polar angles of the measurement direction Xm are specified in a local spherical coordinate system along the reference radial line. In principle, the reference radial line may correspond to any point on the line of sight, but numerical experiments have shown that the selection of the surface point yields the most accurate results. The nadir radiance can then be expressed as Z r etext ðrTOA Þ Iðr s ; X Þjm jrðXmr ; X Þ dO p 2p Z Z rTOA 1 sscat ðrÞ þ Pðr; Xmr ; X0 ÞIðr; X0 Þ dO0 4p mr rs 4p A
IðrA ; Xm Þ ¼
r
r
e½text ðrTOA Þtext ðrÞ dr þ Iss ðrA ; Xm Þ,
(12)
where Iss ðrA ; Xm Þ ¼ Q s ðr s ; Xm Þetext ðjrA rs jÞ þ
Z
J ss ðr; Xm Þetext ðjrA rjÞ ds jrA rs j
is the single scattering radiance and
trext ðrÞ ¼
Z
1
mr
r
sext ðrÞ dr
rs
is the slant extinction optical depth along the reference direction. The above representation shows that the single scattering radiance is calculated for a full spherical geometry, while the multiple scattering radiance is computed for a plane-parallel geometry. Defining the adjoint source function for the nadir radiance by the relation Q y ðr; XÞ ¼
1 sscat ðrÞ A r r r Pðr; Xmr ; XÞe½text ðrTOA Þtext ðrÞ þ dðr r s ÞHðmÞjmjrðXmr ; XÞetext ðrTOA Þ , 4p p
mr
the integral representation Im ¼ hQ y ; Ii þ Iss ,
(13)
with Im ¼ IðrA ; Xm Þ, holds true. It should be pointed out that the equivalence between the representations (13) and (12) follows directly from the definition of the scalar product hQ y ; Ii ¼
Z
Z
rTOA
Q y ðr; XÞIðr; XÞ dO dr,
4p
rs
and the sifting property of the Dirac delta function, which yields A
p
r
etext ðrTOA Þ
Z
rTOA
dðr rs Þ dr rs
¼
A
p
r
etext ðrTOA Þ
Z
Iðr; XÞHðmÞjmjrðXmr ; XÞ dO 4p
Z
Iðr s ; X Þjm jrðXmr ; X Þ dO . 2p
Let us now consider a limb viewing geometry with the line of sight bounded by the near side point A and the far side point A0 , and let the radial coordinate of the tangent point of the line of sight T be r T . The radiance measured by a limb instrument is given by IðrA ; Xm Þ ¼
Z
sscat ðrÞ 4p jrA rA0 j
Z 4p
Pðr; Xm ; X0 ÞIðr; X0 Þ dO0 etext ðjrA rjÞ ds þ Iss ðrA ; Xm Þ,
where Iss ðrA ; Xm Þ ¼
Z
J ss ðr; Xm Þetext ðjrA rjÞ ds. jrA rA0 j
(14)
ARTICLE IN PRESS 468
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
m
A M
r T
r+
M’
rTOA rT
r-
A’
rTOA
Fig. 1. Limb scattering geometry.
The line integral in Eq. (14) can be transformed into a simple integral along the radial coordinate by using the basic result Z Z Z f ðrÞ ds ¼ f ðrþ Þ ds þ f ðr Þ ds jrA rA0 j
jrA rT j rA
Z ¼
rT
jr r j
T A ds ½f ðr Þ þ f ðr Þ ðrÞ dr, dr þ
where rþ and r , with jrþ j ¼ jr j ¼ r, are ‘‘symmetric’’ integration points situated on the near side and on the far side of the line of sight, respectively (Fig. 1), and ds 1 r ðrÞ ¼ ¼ Hðr r T Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . dr mðrÞ 2 r r2 T
As for a nadir viewing geometry, the multiple scattering radiance is computed in a plane-parallel atmosphere and for a reference solar zenith angle. The variation with the solar zenith angle is accounted for in an approximate manner by scaling the reference radiance, that is, Pðr; Xm ; XÞIðr; XÞ cðrÞPðr; Xmr ; XÞIðr; XÞ.
(15)
The scaling factor is computed as sun
cðrÞ ¼
etext ðjrrTOA jÞ , sun etext ðjrr rTOA jÞ
where rr is a vector along the radial line for which it holds true that jrr j ¼ jrj ¼ r. The exponential term expftsun ext ðjr rTOA jÞg reproduces the dependency of the source term on the solar zenith angle and is assumed to act like a multiplicative factor for the radiance fields. Noticing that in practice, the reference calculation corresponds to the tangent point, we obtain Z Z rA ds sscat ðrÞ IðrA ; Xm Þ ¼ Pðr; Xmr ; XÞIðr; XÞ dO0 TðrÞ ðrÞ dr þ Iss ðrA ; Xm Þ, (16) 4p dr rT 4p with TðrÞ ¼ cðrþ Þetext ðjrA r
þ
jÞ
þ cðr Þetext ðjrA r jÞ .
In view of Eq. (16), the adjoint source function for the limb radiance is defined by ds sscat ðrÞ Pðr; Xmr ; XÞTðrÞ ðrÞ, Q y ðr; XÞ ¼ Hðr r T Þ 4p dr in which case, IðrA ; Xm Þ possesses the integral representation (13). We are now in the position to formulate the main result of the adjoint radiative transfer approach: If Iðr; XÞ solves the forward problem
LIðr; XÞ ¼ Q ðr; XÞ,
(17)
Iðr TOA ; X Þ ¼ Iðr s ; Xþ Þ ¼ 0,
(18)
y
and I ðr; XÞ solves the adjoint problem
Ly Iy ðr; XÞ ¼ Q y ðr; XÞ,
(19)
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Iy ðr TOA ; Xþ Þ ¼ Iy ðr s ; X Þ ¼ 0,
469
(20)
then the radiance measured by a satellite instrument can be expressed as ð13Þ
Im ¼ hQ y ; Ii þ Iss ð19Þ
¼ hLy Iy ; Ii þ Iss
ð10Þ;ð18Þ;ð20Þ
¼
hIy ; LIi þ Iss
ð17Þ
¼ hIy ; Q i þ Iss .
(21)
The scalar product representation (21) reveals that the information about the measurement geometry is encapsulated in the adjoint source function Q y (and implicitly in the adjoint radiation field Iy ), while the information about the source conditions are contained in the forward source function Q . 3. Solution of the adjoint radiative transfer equation The forward and the adjoint radiative transfer equations (17) and (19) are related to each other. Replacing X by X in the expression of the adjoint operator Ly, we rewrite the adjoint radiative transfer equation as
Ly ðXÞIy ðr; XÞ ¼ Q y ðr; XÞ. y Defining the conjugate adjoint field ^I by the relation
^Iy ðr; XÞ ¼ Iy ðr; XÞ
(22) 0
0
and using the symmetry property of the phase function Pðr; X; X Þ ¼ Pðr; X ; XÞ, we find the identity [6] y
Ly ðXÞIy ðr; XÞ ¼ L^I ðr; XÞ.
(23)
The solution of the adjoint radiative transfer equation can be found by using the same solution method as for the forward problem with an appropriate source function and boundary conditions. Performing an azimuthal expansion of the radiance field we obtain the boundary-value problems:
m
y Z y y d^Im sscat ðrÞ 1 ðr; mÞ ¼ sext ðrÞ^Im ðr; mÞ þ pm ðr; m; m0 Þ^Im ðr; m0 Þ dm0 2 dr 1 sscat ðrÞ r r þ ð2 dm0 ÞF r pm ðr; m; mr Þe½text ðrTOA Þtext ðrÞ , 4p
^Iy ðr TOA ; m Þ ¼ 0, m ^Iy ðr s ; mþ Þ ¼ 2A m
Z
0
1
^Iy ðr s ; m Þjm jr ðmþ ; m Þ dm þ ð2 dm0 ÞF r A m r ðmþ ; m Þetrext ðrTOA Þ , r m r m m
p
(24)
for a nadir viewing geometry, and
m
y Z y y dI^m sscat ðrÞ 1 ðr; mÞ ¼ sext ðrÞ^Im ðr; mÞ þ pm ðr; m; m0 Þ^Im ðr; m0 Þ dm0 2 dr 1 ds sscat ðrÞ þ ð2 dm0 ÞHðr r T Þ pm ðr; m; mr ÞTðrÞ ðrÞ, 4p dr
^Iy ðr TOA ; m Þ ¼ 0, m ^Iy ðr s ; mþ Þ ¼ 2A m
Z
0
1
^Iy ðr s ; m Þjm jr ðmþ ; m Þ dm m m
(25)
for a limb viewing geometry. In Eqs. (24) and (25), pm and rm are the azimuthal expansion coefficients of the phase function and the bi-directional reflection function, respectively, and F r ¼ 1=mr . It should be pointed out that the definition of the y conjugate adjoint field (22) implies that Iym ðr; mÞ ¼ ^Im ðr; mÞ. The adjoint radiative transfer problems (24) and (25) are solved by using the discrete ordinate method with matrix exponential [13]. The peculiarities of this solution method for the present application will be sketched below. The atmosphere is discretized in N levels: r 1 4r 2 4 4r N with the convention r 1 ¼ r TOA and r N ¼ r s . A layer j, bounded above by the level r j and below by the level r jþ1 , has the geometrical thickness Dr j ¼ r j r jþ1 . The extinction and the scattering coefficients as well as the phase function coefficients are assumed to be constant within each layer; their average values are s¯ ext;j , s¯ scat;j and p¯ m;j , respectively. The angular variation of the phase function and radiance is discretized in M
ARTICLE IN PRESS 470
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
discrete ordinate directions per hemisphere, and the set of Gauss–Legendre quadrature points and weights in the interval ð0; 1Þ is denoted by fmk ; wk gk¼1;M . The matrix exponential formalism operates with the concept of layer equation
A1m;j im;j þ A2m;j im;jþ1 ¼ bm;j ,
(26)
which relates the level values of the adjoint radiance field, 2 1 3 im;j im;j ¼ 4 2 5, im;j 1;2
y
with ½im;j k ¼ ^Im ðr j ; mk Þ and k ¼ 1; . . . ; M. The layer matrices can be expressed in compact form as
A1m;j ¼ D1m;j V1 m;j , A2m;j ¼ D2m;j V1 m;j , where the columns of the matrix Vm;j are the eigenvectors of the layer matrix 2 11 3 Am;j A12 m;j 4 5, Am;j ¼ A12 A11 m;j m;j with entries 1 p¯ ðm ; m Þ 2s¯ ext;j dkl , ½w s¯ 2mk l scat;j m;j k l 1 ½A12 p¯ ðm ; ml Þ, w s¯ m;j kl ¼ 2mk l scat;j m;j k
½A11 m;j kl ¼
and the diagonal scaling matrices D1m;j and D2m;j are given by D1m;j ¼ diag½a0 ðlk Dr j Þ; 1, D2m;j ¼ diag½1; a0 ðlk Dr j Þ, with lk , k ¼ 1; . . . ; M, being the positive eigenvalues of the matrix Am;j and a0 ðxÞ ¼ ex . The expression of the layer vector depends on the viewing geometry. For nadir measurements, bm;j is computed as
bm;j ¼ Lm;j V1 m;j ðDr j bm;j Þ, where the components of the source vector 2 1 3 bm;j bm;j ¼ 4 2 5 bm;j are given by 1;2
½bm;j k ¼
1
mk
ð2 dm0 Þ
Fr p¯ ðmk ; mr Þ. s¯ 4p scat;j m;j
The diagonal scaling matrix Lm;j is expressed as
Lm;j ¼ diag½b1 ðlk Dr j Þ; b2 ðlk Dr j Þ,
(27)
with the interpolation functions b1 and b2 being defined by r
b1 ðxÞ ¼ b2 ðxÞ ¼
r
eðtext;j þxÞ etext;jþ1 , trext;jþ1 trext;j x ðtrext;jþ1 þxÞ
e
trext;j
e
trext;j trext;jþ1 x
.
(28)
The layer equation has been derived by assuming that the extinction optical depth varies linearly within the layers; the level values of the slant extinction optical depth being given by
trext;j ¼
j1 1X
mr u¼1
s¯ ext;u Dru ;
j ¼ 1; . . . ; N.
For a limb viewing geometry, we consider a division of the line of sight in an even number N point LOS of points; the number of point point ¯ limb paths along the line of sight is then Npath LOS ¼ N LOS 1. The index j ¼ N LOS =2 þ 1 corresponds to the tangent point of radial coordinate r T , and for simplicity, we assume that ¯j is also the level index for which it holds true that r ¯j ¼ r T . The
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
471
¯ number of effective layers traversed by the line of sight is Neff lay ¼ j 1, and it is apparent that bm;j vanishes for all layers j4Neff lay . To derive the expression of the layer vector we introduce the following quantities:
(1) the scaling factors sun
cj ¼
etext ðjrj rTOA jÞ ; sun etext ðjrr;j rTOA jÞ
j ¼ 1; . . . ; Npoint LOS ,
where jrj j ¼ jrr;j j ¼ r levðjÞ , and the integer function levðjÞ maps the point index j, into the level index ( j; 1pjp¯j; levðjÞ ¼ ¯ ¯ 2j j; j þ 1pjp2¯j 1; (2) the slant optical depths
trext;1 ¼ 0;
trext;j ¼
j1 X
s¯ ext;layðuÞ Du ; j ¼ 2; . . . ; Npoint LOS ,
u¼1
where Dj are the limb path lengths, and the integer function layðuÞ map the path index u, into the layer index ( u; 1pup¯j 1; layðuÞ ¼ 2¯j u 1; ¯jpup2¯j 2; (3) the average Jacobian values
¯ ds dr
¼ j
Dj ; DrlayðjÞ
j ¼ 1; . . . ; Npath LOS .
The layer vector entering in Eq. (26) can then be expressed as þ 1 1 bm;j ¼ Lþ m;j V m;j ðDr j bm;j Þ þ Lm;j Vm;j ðDr j bm;j Þ, þ
where the components of the source vectors bm;j and bm;j , are given by þ1;2
½bm;j k ¼
1
mk
ð2 dm0 Þ
¯ 1 ds s¯ scat;j p¯ m;j ðmk ; mr Þ¯cj 4p dr j
ð2 dm0 Þ
¯ 1 ds s¯ scat;j p¯ m;j ðmk ; mr Þ¯c2¯jj1 , 4p dr 2¯jj1
and 1;2
½bm;j k ¼
1
mk
respectively. The diagonal matrix Lþ m;j is as in Eqs. (27) and (28), while the interpolation functions b1 and b2 , which enter in r r the expression of the diagonal matrix L m;j , are expressed in terms of the slant optical depths text;2¯jj and text;2¯jj1 in place r r of text;j and text;jþ1 , respectively. The layer equation (26) together with the boundary conditions at the top and the bottom of the atmosphere are assembled into the global matrix of the entire atmosphere, and the solution of the resulting system of equations yields the level values of the adjoint radiance field. An essential aspect of the adjoint approach is the calculation of the scalar product in Eq. (21). The integration with respect to the radial coordinate can be analytically performed (see Appendix B) and this result substantially increases the efficiency of the method. 4. Numerical simulations The purpose of our numerical simulations is to validate the scalar product representation (21), since this equation is the basic relation for computing the weighting functions by the adjoint approach. The radiance measured by a satellite instrument is computed by solving the forward and the adjoint problems in the framework of the discrete ordinate method with matrix exponential. The calculations are performed in the spectral domain between 320 and 350 nm, and for a vertically inhomogeneous atmosphere taking into account Rayleigh and aerosol scattering, as well as the absorption by O3, NO2 and aerosol. The O3 and NO2 profiles have been taken from the model of McLinden et al. [14], while the LOWTRAN/MODTRAN database has been used for aerosol parametrization [15]. The atmosphere is discretized with a step of 1 km between 0 and 25 km, a step of 2.5 km between 25 and 50 km and a step of 5 km between 50 and 100 km. The surface albedo is set to 0.3, and the solar flux is 1.
ARTICLE IN PRESS 472
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Radiance [relative units]
2e−01
1 2 3 4 5 6 7 8
1e−01
0e+00
1 2 3 4 5 6 7 8
320
340
330
350
Relative Difference
2e−03
8 1e−03 8
7 4 3 5 2 1 6
12 3 4 5 6
7
0e+00 320
330
340
350
Wavelength [nm] Fig. 2. (a) Nadir radiances for different values of the solar zenith angle. (b) Relative errors between the solutions computed by using the forward and the adjoint method. The curves correspond to the following values of the solar zenith angle: 1–10 , 2–20 , 3–30 , 4–40 , 5–50 , 6–60 , 7–70 and 8–80 . The local viewing zenith angle is 10 and all angles refer to the top of the atmosphere.
The results shown in Fig. 2 correspond to a nadir viewing geometry. In Fig. 2a we plot the radiances computed by the forward method for different values of the solar zenith angle. The relative errors between the forward and the adjoint methods are smaller than 0.13% (Fig. 2b), while the computational costs are similar. In Fig. 3 we illustrate the simulation results for a limb viewing geometry. For the forward method, the radiative transfer equation is solved for five values of the solar zenith angle (five right-hand side vectors of the matrix equation), while, in view of the approximation (15), a single calculation is performed for the adjoint method. Due to the approximation employed, the relative errors are of about 0.8%. The computer time needed to simulate the radiances at 121 wavelengths and eight discrete ordinate directions per hemisphere was 24 s for the forward method and 14 s for the adjoint method (on an Intel Xeon processor with 2.8 GHz). The main problem of the adjoint method is that for large optical depths along the line of sight when the tangent point is close to the terminator, the approximation (15) deteriorates and the relative errors may reach 3% (Fig. 4). However, in atmospheric remote sensing, due to the general nature of the Gauss–Newton iterative method commonly used to solve the nonlinear inverse problem, accuracies of the weighting functions up to a few percent do not deteriorate the convergence rate of the solution, and the application of the adjoint approach is not critical.
5. Conclusions The adjoint radiative transfer is formulated for a pseudo-spherical atmosphere and general viewing geometries. The single scattering radiance is computed in a spherical atmosphere, while for the multiple scattering radiance we formulate an adjoint radiative transfer equation in a plane-parallel atmosphere starting from the integral representation of the radiance measured by a satellite instrument. For limb scattering geometries, we solve the one-dimensional radiative transfer problem for a reference solar zenith angle, and scale the calculated radiances by a correction factor depending on the solar optical depth. The solution of the adjoint radiative transfer equation is obtained by using the discrete ordinate method with matrix exponential. The accuracy of the adjoint model is satisfactory for a nadir viewing geometry and still acceptable for a limb viewing geometry. The adjoint approach together with a linearized forward model are implemented in a software package that we developed at the German Aerospace Centre for solving inverse problems in atmospheric remote sensing. A comparison of these two linearization methods with respect to their efficiency and accuracy will be the subject of a future paper.
ARTICLE IN PRESS
Radiance [relative units]
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
473
2e−01
96 5 78
1e−01
4 3
9 8 7 6 5 4 3
0e+00 12 320
12
330
340
350
340
350
Relative Difference
8e−03 9 8 7 4 6e−03 56
9 8 6 45 7
4e−03 3
3 12
2e−03 2 1
0e−00 320
330 Wavelength [nm]
Fig. 3. (a) Limb radiances (b) Relative errors between the solutions computed by using the forward and the adjoint method. The curves correspond to the following values of the tangent height: 1–50, 2–40, 3–30, 4–20; 5–16; 6–13; 7–10; 8–6; 9–3 km. The solar zenith angle is 45 and the relative azimuthal angle is 45 .
Radiance [relative units]
2e−01 98765
1e−01 4
0e+00
3
5 7 6 9 8 4 3 2 1
320
21
330
340
350
340
350
Relative Difference
3e−02 9 8
2e−02
9 8
7 6 5
7
1e−02
4
6
3 1 2
5 12 3 4
0e+00 320
330 Wavelength [nm]
Fig. 4. Same as in Fig. 3 but for a solar zenith angle of 80 and a relative azimuthal angle of 0 .
ARTICLE IN PRESS 474
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
Acknowledgement The authors are grateful to Eugene Ustinov for illuminating discussions regarding the modeling of limb radiances in the framework of adjoint radiative transfer. Appendix A In this appendix we prove that a diffuse radiation field solving the inhomogeneous operator equation (8) with the boundary conditions (9) also solves the radiative transfer equation (1) with the boundary conditions (2) and (3). The proof relies on the equivalence between the boundary-value problem (1)–(3) and the integral representations of the radiative transfer equation for upwelling and downwelling radiation: If the diffuse radiation field Iðr; XÞ satisfies the boundary-value problem (1)–(3), then Iðr; XÞ solves the integral equation Z rTOA 1 0 Iðr; X Þ ¼ Jðr 0 ; X Þe½text ðr Þtext ðrÞ dr 0 (29) jm j r for X ¼ X, with
text ðrÞ ¼
Z
1 jm
j
r
sext ðrÞ dr rs
and þ
Iðr; Xþ Þ ¼ Iðr s ; Xþ Þetext ðrÞ þ
1
m
þ
Z
r
þ
þ
Jðr 0 ; Xþ Þe½text ðrÞtext ðr Þ dr 0 0
(30)
rs
for X ¼ Xþ, with
tþext ðrÞ ¼
1
m
þ
Z
r
sext ðrÞ dr. rs
The converse results is also true. Let us now assume that the diffuse radiation field Iðr; XÞ satisfies the boundary-value problem (8)–(9). For downwelling radiances, there hold X ¼ X and m ¼ m o0. As a result we deduce that HðmÞ ¼ 0, and the differential equation (8) can be written explicitly as
m
dI ðr; X Þ ¼ sext ðrÞIðr; X Þ þ Jðr; X Þ. dr
(31)
Integrating the inhomogeneous differential equation (31) with the boundary condition Iðr TOA ; X Þ ¼ 0; yields the integral form of the radiative transfer equation for the downwelling radiation (29). For upwelling radiances, the conditions X ¼ Xþ and m ¼ mþ 40 imply that HðmÞ ¼ 1, and the differential equation (8) read as Z dI A mþ ðr; Xþ Þ ¼ sext ðrÞIðr; Xþ Þ þ Jðr; Xþ Þ þ dðr rs Þmþ Iðr; X Þjm jrðXþ ; X Þ dO dr p 2p þ dðr r s Þmþ Q s ðr; Xþ Þ.
(32) þ
Integrating the inhomogeneous differential equation (32) with the boundary condition Iðr s ; X Þ ¼ 0, and using the results Z r Z þ þ 1 A 0 dðr 0 rs Þmþ Iðr 0 ; X Þjm jrðXþ ; X Þ dO e½text ðrÞtext ðr Þ dr 0 þ m rs p 2p Z þ A ¼ Iðr s ; X Þjm jrðXþ ; X Þ dO etext ðrÞ
p
2p
and 1
mþ
Z
r rs
þ
þ
þ
½dðr 0 r s Þmþ Q s ðr 0 ; Xþ Þe½text ðrÞtext ðr Þ dr 0 ¼ Q s ðr s ; Xþ Þetext ðrÞ , 0
we find the integral form of the radiative transfer equation for the upwelling radiation (30). Appendix B In this appendix we sketch the calculation of the scalar product in Eq. (21) for a nadir viewing geometry; the treatment for a limb viewing geometry is similarly. Considering the multiple scattering radiance Ims ¼ hIy ; Q i, and assuming the
ARTICLE IN PRESS A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
475
azimuthal expansions 2M1 X Ims ðXr Þ ¼ Ims;m ðmr Þ cos mjr , m¼0
Iy ðr; XÞ ¼
2M1 X
Iym ðr; mÞ cos½mðj jr Þ,
m¼0
Q ðr; XÞ ¼
2M1 X
Q m ðr; mÞ cos mj,
m¼0
we find that Ims;m ¼ ð1 þ dm0 Þp
Z
r TOA rs
Z
1
Q m ðr; mÞIym ðr; mÞ dm dr.
1
Using the expression of the mth Fourier harmonics of the source function, sscat ðrÞ sun Q m ðr; mÞ ¼ ð2 dm0 ÞF sun pm ðr; m; msun Þetext ðjrrTOA jÞ 4p A sun þ ð2 dm0 ÞF sun dðr r s ÞHðmÞmjmsun jrm ðm; msun Þetext ðjrs rTOA jÞ ,
p
we obtain 1 sun 1 NX 2 F sun sTm;j Im;j þ 2AF sun jmsun jetext;N rTm im;N , 2 j¼1
Ims;m ¼
where the components of the scattering vector sm;j are given by ¼ wl s¯ scat;j p¯ m;j ð ml ; msun Þ, ½s1;2 m;j l the reflection vector rm has the entries ½rm l ¼ wl ml rm ðml ; msun Þ, and accordingly to the average secant approximation, Im;j is the integral Z Dr j ½ðr=Dr j Þtsun þð1r=Drj Þtsun ext;j ext;jþ1 i Im;j ¼ e (33) m;j ðrÞ dr. 0
To compute Im;j we use the analytical representation of the internal radiance derived in [13], that is, 2 3 " # 1 diag½elk ðDrj rÞ xm;j diag½B1 ðlk ðDr j rÞÞ 0 5 þ Vm;j V1 Dr b , im ðrÞ ¼ Vm;j 4 2 0 diag½B2 ðlk rÞ m;j j m;j diag½elk r xm;jþ1 1
2
1
2
1 T with xm;j ¼ ½xm;j ; xm;j T ¼ V1 m;j im;j , xm;jþ1 ¼ ½xm;jþ1 ; xm;jþ1 ¼ V m;j im;jþ1 , and ðtrext;j þxÞ
B1 ðxÞ ¼
e
½ðr=Dr j Þtrext;j þð1r=Dr j Þtrext;jþ1
e
r
B2 ðxÞ ¼
,
trext;jþ1 trext;j lk Drj r
r
eðtext;jþ1 þxÞ e½ðr=Drj Þtext;j þð1r=Drj Þtext;jþ1 . trext;j trext;jþ1 lk Drj
The result of the integration is 2 3 " 1 ½diagða1 ðlk Dr j ÞÞxm;j ½diagðB¯ 1 ðlk Dr j ÞÞ 5 þ Vm;j Im;j ¼ Vm;j 4 2 0 ½diagða2 ðlk Dr j ÞÞxm;jþ1
0 ½diagðB¯ 2 ðlk Dr j ÞÞ
# V1 m;j Dr j bm;j ,
where sun a1 ðxÞ ¼ Dr j I0 ðtsun ext;j ; text;jþ1 þ xÞ, sun a2 ðxÞ ¼ Dr j I0 ðtsun ext;j þ x; text;jþ1 Þ,
B¯ 1 ðxÞ ¼ B¯ 2 ðxÞ ¼
Drj
trext;jþ1 trext;j lk Dr j Drj
trext;j trext;jþ1 lk Dr j
trext;j
½e
r
sun sun r sun r ½etext;jþ1 I0 ðtsun ext;j þ x; text;jþ1 Þ I0 ðtext;j þ text;j ; text;jþ1 þ text;jþ1 Þ,
and I0 is the integral Z 1 I0 ða; bÞ ¼ e½axþbð1xÞ dx ¼ 0
sun sun r sun r I0 ðtsun ext;j ; text;jþ1 þ xÞ I0 ðtext;j þ text;j ; text;jþ1 þ text;jþ1 Þ,
1 ðea eb Þ. ab
References [1] Marchuk GI. Equation for the value of information from weather satellites and formulation of inverse problems. Cosmic Res 1964;2: 394–409.
ARTICLE IN PRESS 476
A. Doicu, T. Trautmann / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 464–476
[2] Box MA, Gerstl SAW, Simmer C. Application of the adjoint formulation to the calculation of atmospheric radiative effects. Beitr Phys Atmos 1988;61:303–11. [3] Ustinov EA. Inverse problem of photometric observation of solar radiation reflected by an optically dense planetary atmosphere: mathematical methods and weighting functions of linearized inverse problem. Cosmic Res 1991;29:519–32. [4] Ustinov EA. Adjoint sensitivity analysis of radiative transfer equation: temperature and gas mixing ratio weighting functions for remote sensing of scattering atmospheres in thermal IR. JQSRT 2001;68:195–211. [5] Ustinov EA. Atmospheric weighting functions and surface partial derivatives for remote sensing of scattering planetary atmospheres in thermal spectral region: general adjoint approach. JQSRT 2005;92:351–71. [6] Rozanov VV, Rozanov AV. Generalized form of the direct and adjoint radiative transfer equations. JQSRT 2007;104:155–70. [7] Landgraf J, Hasekamp OP, Box MA, Trautmann T. A linearized radiative transfer model for ozone profile retrieval using the analytical forward-adjoint perturbation theory approach. J Geophys Res 2001;106(D21):27291–305. [8] Walter H, Landgraf J, Hasekamp OP. Linearization of a pseudo-spherical vector radiative transfer model. JQSRT 2004;85:251–83. [9] Ustinov EA. Adjoint approach to evaluation of weighting functions for remote sensing of scattering planetary atmospheres in thermal spectral region with limb-viewing geometry. EGU General Assembly 2008, Vienna, Geophysical Research Abstracts, vol. 10, EGU2008-A-10664, 2008, SRef-ID: 16077962/gra/EGU2008-A-10664. [10] Rozanov VV, Kurosu T, Burrows JP. Retrieval of atmospheric constituents in the UV-visible: A new quasi-analytical approach for the calculation of the weighting functions. JQSRT 1998;60:277–99. [11] Rozanov VV, Kokhanovsky AA. The average number of photon scattering events in vertically inhomogeneous atmosphere. JQSRT 2005;96:11–33. [12] Qin Y, Jupp DLB, Box MA. Extension of the discrete-ordinate algorithm and efficient radiative transfer calculation. JQSRT 2002;74:767–81. [13] Doicu A, Trautmann T. Discrete ordinate method with matrix exponential for a pseudo-spherical atmosphere. Scalar case. JQSRT 2009;110:146–58. [14] McLinden CA, McConnell JC, Griffioen E, McElroy CT. A vector radiative-transfer model for the Odin/OSIRIS project. Can J Phys 2002;80:375–93. [15] Kneizys F, Abreu L, Anderson G, Chetwynd J, Shettle E, Berk A, Bernstein L, Robertson D, Acharya P, Rothman L, Selby J, Gallery W, Clough S. The MODTRAN 2/3 report and LOWTRAN 7 model. Technical Report, contract F19628-91-C-0132 with Ontar Corporation, Phillips Laboratory, Hanscom AFB; 1996.