Towards linearization of atmospheric radiative transfer in spherical geometry

ARTICLE IN PRESS

Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200 www.elsevier.com/locate/jqsrt

Towards linearization of atmospheric radiative transfer in spherical geometry Holger H. Walter, Jochen Landgraf National Institute for Space Research (SRON), Earth Oriented Science Division, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands Received 30 January 2004; accepted 31 August 2004

Abstract We present a general approach for the linearization of radiative transfer in a spherical planetary atmosphere. The approach is based on the forward-adjoint perturbation theory. In the first part we develop the theoretical background for a linearization of radiative transfer in spherical geometry. Using an operator formulation of radiative transfer allows one to derive the linearization principles in a universally valid notation. The application of the derived principles is demonstrated for a radiative transfer problem in simplified spherical geometry in the second part of this paper. Here, we calculate the derivatives of the radiance at the top of the atmosphere with respect to the absorption properties of a trace gas species in the case of a nadir-viewing satellite instrument. r 2004 Elsevier Ltd. All rights reserved. Keywords: Radiative transfer; Spherical geometry; Linearization; Adjoint formulation; Perturbation theory

1. Introduction Current and future ultraviolet (UV) and visible (VIS) satellite instruments, which observe the Earth’s atmosphere in limb viewing geometry, provide information about the distribution of atmospheric trace constituents with high vertical resolution. For example the Scanning Imaging Corresponding author. Tel.: +31 30 2535913; fax: +31 30 2540860.

E-mail address: [email protected] (HH. Walter). 0022-4073/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2004.08.043

ARTICLE IN PRESS 176

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) [1], launched in 2002 on board of ESA’s ENVISAT satellite, measures backscattered limb radiances in the wavelength range from 240 to 2380 nm. In limb viewing geometry, SCIAMACHY vertically scans the atmosphere from ground level up to an altitude of approximately 100 km in steps of 3 km at tangent height. This enables the retrieval of a variety of atmospheric trace constituents in high vertical resolution. Furthermore, SCIAMACHY offers the unique possibility to observe an air volume in nadir as well as in limb viewing geometry. It is expected that the combination of both viewing geometries will facilitate a higher accuracy in the retrieval of atmospheric trace gas profiles in the troposphere and stratosphere [2]. In 2001 the Canadian Optical Spectrograph and Infrared Imager System (OSIRIS) [3] was launched on board the Odin satellite. The instrument observes the backscattered radiances in limb viewing geometry and comprises a UV/ VIS spectrograph covering the wavelength range from 280 to 800 nm and three 30 nm wide infrared (IR) channels at 1260, 1270 and 1280 nm. The vertical distribution of e.g. O3 ; NO2 ; OClO, BrO and aerosols can be retrieved from OSIRIS limb measurements. Also, the Ozone Mapping and Profiler Suite (OMPS), to be launched on NPOESS in 2010 by the United States of America, is dedicated to monitoring the global ozone distribution from space. This instrument consists of a nadir sensor with a wide field-of-view telescope and a limb sensor. The nadir spectrometer observes the backscattered radiance between 300 and 380 nm, whereas the limb sensor measures the along-track limb-scattered solar radiance in the spectral range from 290 to 1000 nm. In order to retrieve an atmospheric parameter x from such limb observations with standard least squares methods [4–7], one has to simulate both the radiance I TOA at the top of the model atmosphere (TOA) and its derivative with respect to the parameter to be retrieved. This corresponds to a linearization of the radiative transfer problem, given by a Taylor expansion around a first-guess parameter x0 : I TOA ðxÞ ¼ I TOA ðx0 Þ þ

qI TOA ðx0 ÞDx þ OðDx2 Þ: qx

(1)

Here, OðDx2 Þ represents higher order terms in Dx ¼ x  x0 : For the simulation of limb radiance measurements in the UV and VIS spectral range, radiative transfer models, which take multiple scattering of radiation in a spherical model atmosphere into account, are required. In general, such spherical radiative transfer models are complex and computationally highly expensive. The Monte Carlo technique [8] provides a straightforward approach to solving this type of radiative transfer problem (see e.g. [9–11]). However, due to its statistical nature, Monte Carlo models require a considerable computational effort in order to achieve a reasonable accuracy. Therefore alternative techniques have been developed in order to solve the spherical radiative transfer equation in a more efficient manner. Sen and Wilson [12] give a general overview about different solution techniques in curved media. Here, applications such as radiative transfer in stellar atmospheres are discussed. For the simulation of radiative transfer in a spherical planetary atmosphere, Sobolev and Minin [13,14] suggested an efficient solution approach, in which the first-order scattering is calculated exactly, whereas the higher-order scattering is approximated by the use of the first two terms of a Legendre expansion of the scattering phase function only. Smokty [15] applied this method to determine the source function

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

177

and brightness of a spherical planetary atmosphere. Furthermore, several radiative transfer approaches have been suggested [16–19] for the simulation of the azimuthally-averaged radiation field in a spherical atmosphere. However, remote sensing in limb-viewing geometry requires the simulation of the azimuthal dependence of the radiation field. This angular dependence is fully taken into account by the model of Herman et al. [20], which employs a Gauss–Seidel iteration scheme to solve the integral form of the radiative transfer equation. Here, the assumption of a constant ratio of single scattering to multiple scattering along a conical boundary facilitates the efficient calculation of the multiple-scattering source function in spherical geometry. However, due to this approximation the radiation field can only be calculated in cases in which the single-scattering terms are not zero. Rozanov et al. [21] presented recently an alternative radiative transfer model, which is based on a Picard iteration scheme. A radiation field, which is calculated in pseudo-spherical approximation, serves as an initial guess in order to integrate the radiative transfer equation along discrete lines in spherical geometry in subsequent iterations. Minor attention has been given to the linearization of spherical radiative transfer. Kaiser and Burrows [22] have presented a linearized model, taking two orders of scattering into account. However, in order to linearize a multiple-scattering radiative transfer problem in spherical geometry, the finite difference method has been employed so far. This means that two independent radiative transfer calculations are needed to approximate the derivative of I TOA in Eq. (1): one for an unperturbed atmosphere and one for an atmosphere for which the parameter x0 has been perturbed by a small amount. For example, for the retrieval of a vertical trace gas distribution such a method obviously requires lots of repetitive computations. Each perturbation of the trace gas concentration at a given altitude requires a separate radiative transfer calculation. This method combined with the costly spherical radiative transfer models results in an enormous computational effort especially for the application in an iterative retrieval approach. Hence, there is a clear need for a (quasi-)analytical linearization of full spherical radiative transfer. A general linearization approach is provided by the forward-adjoint perturbation theory, known from neutron transport theory (see e.g. [23]). The application of perturbation theory in the context of atmospheric inversion theory was first proposed by Marchuk [24], who applied it to plane-parallel atmospheric radiative transfer. This technique allows the linearization of any modeled observation of the radiation field with respect to absorption and scattering properties of the atmosphere. Miscellaneous applications of the forward-adjoint perturbation theory in planeparallel geometry are discussed e.g. by Box et al. [25], Ustinov [26–28], Rozanov et al. [29], Landgraf et al. [30] and Hasekamp and Landgraf [31]. Due to the fact that the perturbation theory approach can be formulated for an arbitrary geometry [23], it can be applied also to the full spherical radiative transfer problem, which is the subject of this paper. In Section 2, the forward and the adjoint formulations of radiative transfer in spherical geometry are introduced. Both formulations of radiative transfer are needed in Section 3 to derive an analytical linearization of the spherical radiative transfer problem. Any spherical radiative transfer model, which simulates the internal forward and adjoint radiation field, can be linearized by means of the approach presented here. Section 4 discusses symmetries which can be used to simplify the calculation of the forward and adjoint radiation field. In Section 5 the linearization with respect to the amount of trace gases at different altitudes is demonstrated for a

ARTICLE IN PRESS 178

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

simplified spherical radiative transfer problem, which assumes an azimuthally-independent radiation field.

2. Radiative transfer in spherical geometry 2.1. The radiative transfer problem in its forward formulation The radiance and the polarization of radiation can be described by an intensity vector I; which has the Stokes parameters I, Q, U and V as its components I ¼ ½I; Q; U; V T ;

(2)

where T indicates the transposed vector [32]. A common approach in radiative transfer is to neglect the polarization components Q, U and V, and to characterize the radiation field by its radiance component I only. Although this biases the simulation of radiative transfer [32,33], it greatly simplifies the numerical simulation without affecting the basic structure of the radiative transfer equation. Thus, in order to demonstrate the linearization principles of spherical radiative transfer with the forward-adjoint perturbation theory approach, we employ the scalar representation of radiative transfer. The formulation can be extended to the full vector radiative transfer representation in a straightforward manner as already demonstrated for a plane-parallel model atmosphere [31]. The time-independent scalar radiative transfer equation is a detailed balance equation, which describes the gain or loss of radiation along some path through a participating medium. It is given by [32] X  =Iðr; XÞ ¼ be ðrÞI þ Sðr; XÞ;

(3)

where I is the spectral radiance, be is an extinction coefficient and S indicates the source function describing emission and scattering of light. The scalar product X  = denotes the derivative in direction X; where = is the nabla operator. The radiation field depends on the position vector r and on the direction X of the radiation at the position r: The domain of both variables r and X constitutes the so-called phase space on which the radiative transfer problem is defined. For radiative transfer in a planetary atmosphere it is common to represent the position vector r in spherical coordinates, such as the radius r, the cosine of the global zenith angle cos C ¼ Z and the global azimuth angle F: The vector X ¼ ðm; jÞ is specified by the cosine of the local zenith angle y; m ¼ cos y; and the local azimuth angle j: Due to the specific choice of the coordinate system, which is illustrated in Fig. 1, a constant direction is represented by local coordinates y and j; which depend on the position r (compare e.g. [21]). In the remainder of this paper, this dependence will not be explicitly expressed. The first term in Eq. (3) describes the change of radiation in a certain direction and is called ‘streaming term’. The second term describes the attenuation of radiation due to extinction processes. Finally, the third term of Eq. (3) comprises the source function S and therefore represents the gain of photons due to scattering from all directions into the direction X and due to thermal emission of the atmosphere. For radiative transfer in the ultraviolet and visible spectral range, thermal emission can be neglected. Therefore, the source

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

179

solar radiation . . .

. . .

Z θ0

z θ x Ψ

r

ϕ

Ω y

Φ

Y

X

Fig. 1. Schematic overview of the geometry used for the calculation of radiative transfer in a spherical atmosphere. In a spherical coordinate system the radiation field is described by the three spatial coordinates r; C; F and by the two directional coordinates y and j: The incident solar radiation illuminates one hemisphere of the planet. Here, the solar zenith angle y0 is defined locally at the top of the model atmosphere and is changing continuously over the hemisphere.

function S is given by Z bs ðrÞ ~ XÞIðr; XÞ ~ dO ~ Pðr; X; Sðr; XÞ ¼ 4p 4p Z A ~ XÞIðr; XÞ ~ dO ~ þ Hðr; X; 4p p þ S0 ðr; XÞ

ð4Þ

~ ¼ dm~ dj: ~ The first term on the right-hand side of Eq. (4) depicts the gain of photons due with dO to scattering, where bs is a scattering coefficient and P is the scattering phase function. The second term represents the reflection of photons at the bottom of the model atmosphere (BOA) with Lambertian surface albedo A. Here, ~ XÞ ¼ dðr  rBOA ÞYðmÞjmjYðmÞj ~ mj ~ Hðr; X;

(5)

denotes a reflection function. d is the Dirac-delta function and Y represents the Heaviside step function. The Earth radius, i.e. the distance from the center of the Earth to the lower boundary, is

ARTICLE IN PRESS 180

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

given by rBOA : The third term of Eq. (4), S0 ; describes external radiation sources. In the case of solar radiation illuminating one hemisphere of the Earth’s atmosphere (see Fig. 1), S 0 is given by S0 ðr; XÞ ¼ mF 0 dðr  rTOA ÞYðm0 ÞdðX  X0 Þ:

(6)

Here, F 0 is the extraterrestrial flux, rTOA is the distance from the center of the Earth to the top of the model atmosphere and X0 ¼ ðm0 ; j0 Þ describes the direction of the solar beam, with m0 ¼ cos y0 : Here, y0 is the local solar zenith angle and j0 the local solar azimuth angle, both defined at the top of the atmosphere. Further, dðX  X0 Þ ¼ dðm  m0 Þ dðj  j0 Þ: The Heaviside step function Yðm0 Þ in Eq. (6) accomplishes the illumination of only one hemisphere of the Earth’s atmosphere. Due to the fact that the solar illumination as well as the surface reflection are included in the source function S, the radiation field I is subject to the free surface boundary conditions of no incoming radiation at the top and bottom of the atmosphere, IðrTOA ; XÞ ¼ 0

for

IðrBOA ; XÞ ¼ 0

for 0omp1:

 1pmo0; ð7Þ

In the context of remote sensing a planetary atmosphere one is generally interested in the simulation of the radiance at the top of the model atmosphere in the viewing direction of the instrument IðrTOA ; Xv Þ I TOA : Here, Xv ¼ ðmv ; jv Þ describes the viewing direction of the satellite instrument, where mv is the cosine of the viewing zenith angle yv and jv is the viewing azimuth angle, both defined at the location of the satellite. I TOA may be interpreted as a special form of the radiative effect E, defined by the inner product [24,25] E ¼ hRjIi;

(8)

where the inner product describes an integration over the full phase space [23], viz Z Z hRjIi ¼ dV dO Rðr; XÞIðr; XÞ: V

(9)

4p

Here, dV represents a differential volume element. In spherical coordinates it is given by dV ¼ r2 dr dZ dF: R is the so-called response function defined by the measurement. For a measurement of I TOA with a satellite, the response function is given by the direction of the instrument’s line-ofsight Rðr; XÞ ¼ dðr  rTOA ÞdðX  Xv Þ

(10)

with dðr  rTOA Þ ¼

1 dðr  rTOA Þ dðZ  ZTOA Þ dðF  FTOA Þ r2

(11)

in the representation of spherical coordinates. This formal notation of the measurement simulation is essential in order to derive the linearization of radiative transfer in a spherical planetary atmosphere employing the forward-adjoint perturbation theory. Here, the solution of a so-called adjoint transfer equation is needed in addition, which is the subject of the following section.

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

181

2.2. The adjoint formulation of radiative transfer The adjoint transfer equation corresponds to a backward formulation of radiative transfer in space and can be derived from the forward radiative transfer equation (3). In order to derive the adjoint transfer equation, we will follow the theoretical approach introduced by Bell and Glasstone [23]. Eq. (3) can be written in operator form [34] as LI ¼ S 0

(12)

with the transport operator   Z   b ðrÞ A s ~ ~ ~ ~ X  = þ be ðrÞ dðX  XÞ  dO Pðr; X; XÞ  Hðr; X; XÞ : L¼ 4p p 4p

(13)

The operator Ly ; adjoint to L; is defined by the requirement that the inner product fulfills [23–25] hI y jLIi ¼ hLy I y jIi:

(14)

y

Here, I denotes the adjoint intensity field. Just as the intensity I satisfies the free surface boundary conditions, Eq. (7), of no incoming radiation at the top and the bottom of the atmosphere, so too the adjoint intensity I y has to satisfy the boundary conditions of no outgoing adjoint intensity at the top and the bottom of the atmosphere, I y ðrTOA ; XÞ ¼ 0 y

I ðrBOA ; XÞ ¼ 0

for 0pmo1; for  1omp0:

In accordance with Eq. (14) the adjoint transport operator may then be defined as   Z   bs ðrÞ A y ~ ~ ~ ~ Pðr; X; XÞ  Hðr; X; XÞ : L ¼ dO X  = þ be ðrÞ dðX  XÞ  4p p 4p

ð15Þ

(16)

In Appendix A it is shown that the transport operator Ly indeed fulfills Eq. (14). The only difference of the adjoint transport operator, Eq. (16), to the forward transport operator, Eq. (13), is a minus sign in the streaming term and an exchange of directions in the scattering phase function P and in the reflection function H. Analogous to Eq. (12) it is possible to define an adjoint transfer equation Ly I y ¼ Sy0 ;

(17)

where Sy0 denotes any suitable adjoint radiation source. The adjoint transfer problem is closely linked to the forward radiative transfer problem via Eq. (14). Taking Eqs. (12) and (17) into account yields hI y jS 0 i ¼ hS y0 jIi:

(18)

A comparison of Eq. (18) with Eq. (8) reveals that a proper choice of the adjoint source offers the possibility to calculate the radiative effect E in two alternative ways [25]. If we choose the response function, Eq. (10), as an adjoint source Sy0 ðr; XÞ ¼ Rðr; XÞ;

(19)

ARTICLE IN PRESS 182

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

then the radiative effect E can be expressed in two ways: E ¼ hRjIi ¼ hI þ jS0 i:

(20)

This shows that the solution of the adjoint transfer equation is complementary to the solution of the forward radiative transfer equation and has a clear physical meaning. The adjoint intensity I þ can be interpreted as the importance [23,35] of photons within the atmosphere for a given measurement. Thus, in our case, where E represents the radiance I TOA at the top of the model atmosphere in viewing direction of the satellite instrument, the importance I þ in Eq. (20) measures the contribution of photons originating from a radiation source S0 to the satellite measurement. The adjoint transfer problem may be transformed into a pseudo-forward problem [23], which allows one to utilize the same radiative transfer model for the solution of the forward radiative transfer equation (12) as well as for the solution of the adjoint transfer equation (17). For this purpose the phase function has to meet the symmetry relation Pðr; X1 ; X2 Þ ¼ Pðr; X2 ; X1 Þ;

(21)

which is satisfied if P is a function of the scattering angle only [36]. Then, with the substitution of the adjoint intensity field I y by the pseudo-forward intensity field x; defined by xðr; XÞ ¼ I y ðr; XÞ; the adjoint formulation, Eq. (17), can be written as Z bs ðrÞ 0 0 0 Pðr; X0 ; X00 Þxðr; X00 Þ dO00 X  =xðr; X Þ þ be ðrÞxðr; X Þ  4p 4p Z A Hðr; X0 ; X00 Þxðr; X00 Þ dO00 ¼ Rðr; X0 Þ  p 4p

(22)

ð23Þ

~ in the integral terms with a change of the variable X into X0 ; and a change of the variable X 00 into X as a new variable of integration. Due to the symmetry relation of the phase function P, Eq. (21), and an analogous relation for the reflection function H, Eq. (23) can be written in pseudo-forward formulation as Lx ¼ S x

(24)

Sx ðr; XÞ ¼ Rðr; XÞ

(25)

with

as a new radiation source [23]. The transport operator in Eq. (24) is identical to the forward transport operator L from Eq. (13). Additionally, the boundary conditions from Eq. (15) transform to boundary conditions analogous to Eq. (7). For further details we refer to Bell and Glasstone [23].1 With the transformation above it is possible to use the same radiative transfer model for the solution of the forward as well as of the pseudo-forward transfer equation, Eqs. (12) and (24), 1

In the adopted pseudo-forward formulation, Eq. (18) corresponds now to the concept of optical reciprocity [37], which shows again the close interconnection of the two transfer problems.

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

183

respectively. However, it should be noted that there is a substantial difference in the definition of the forward solar source, Eq. (4), and the pseudo-forward radiation source, Eq. (25). The forward solar source S 0 illuminates the complete hemisphere of the Earth’s atmosphere whereas the pseudo-forward radiation source Sx is defined at the satellite’s location only and points into the viewing direction of the instrument. The latter case defines a so-called searchlight problem (compare e.g. [38]). In general, three-dimensional radiative transfer calculations are required to solve the searchlight problem. However, in the case of intrinsic symmetries of the forward radiation field, the calculation of the associated pseudo-forward radiation field may be simplified substantially. This subject will be discussed in more detail in Section 4.

3. Perturbation theory in spherical geometry With the operator formulation of the forward and the adjoint transfer problem from Section 2 it is possible to derive the linearization principles appropriate for full spherical radiative transfer. Suppose I 0 is the solution of the forward radiative transfer problem, Eq. (12), and I y0 is the solution of the adjoint transfer problem, Eq. (17), for a given model atmosphere, which is characterized by a parameter x0 : A perturbation of x0 by an amount Dx to x ¼ x0 þ Dx causes a change of the transport operator L by DL ¼ LðxÞ  Lðx0 Þ; but leaves the radiation source S 0 in Eq. (4) and the response function R in Eq. (10) unaffected. This in turn results in a perturbation of the radiative effect, viz EðxÞ ¼ Eðx0 Þ þ DE:

(26)

Here, the perturbation DE can be calculated from the forward and the adjoint intensity field I 0 and I y0 to first order in Dx by DE ¼ hI y0 jDLI 0 i þ OðDx 2 Þ;

(27)

where OðDx2 Þ denotes second and higher order terms. The derivation of Eq. (27) is analogous to that presented by Marchuk [24] and Box [25,39] for a plane-parallel atmosphere. Then, a comparison of Eq. (26) with a corresponding first-order Taylor expansion of EðxÞ around x0 provides an analytical expression for the derivative qE 1 y ðx0 Þ ¼  hI jDLI 0 i: (28) qx Dx 0 The strength of this approach becomes obvious if derivatives of E are needed with respect to several atmospheric parameters xi with i ¼ 1; . . . ; N: The perturbation of each of these parameters causes a different perturbation DLi of the radiative transport operator. However, for the corresponding derivatives qE=qxi in Eq. (28), the forward and the adjoint intensity field stay the same in all cases. Therefore the calculation of all those derivatives requires only the evaluation of the inner product in Eq. (28). Depending on the number of parameters N this results in an efficient linearization approach for spherical radiative transfer compared to a finite difference scheme. The approach presented here can be used to linearize any spherical radiative transfer model that calculates the internal forward and adjoint intensity fields. The linearization can be achieved with

ARTICLE IN PRESS 184

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

respect to absorption and scattering properties of the model atmosphere. For example, if x represents the number density of an absorbing trace gas in an atmospheric volume DV the perturbation in the transport operator DL is given by DL ¼ DxGðrÞsabs ðrÞ

(29)

with ( GðrÞ ¼

1;

r 2 DV ;

0;

reDV

(30)

and sabs ðrÞ is the absorption cross section of the considered trace gas species. Substitution of Eq. (29) in Eq. (28), whilst taking Eq. (22) into account, yields the derivatives Z Z qE dV sabs ðrÞ dO x0 ðr; XÞI 0 ðr; XÞ: (31) ðx0 Þ ¼  qx DV 4p Thus, Eq. (31) provides the linearization of a spherical radiative transfer model with respect to a trace gas density in its most general form, by identifying the radiative effect E with the radiance at the top of the model atmosphere I TOA in Eq. (1). In the derivation of Eq. (28) and subsequently of Eq. (31) we made use of the fact that the radiation source S 0 in Eq. (4) and the response function R in Eq. (10) are defined at the top of the model atmosphere and are not affected by any perturbation Dx inside the atmosphere. Thus, DS0 ¼ 0 and DR ¼ 0: For cases in which DS 0 and DR do not vanish the corresponding derivation of the derivatives may be found in [26,40].

4. Symmetry properties Symmetry properties of a model atmosphere combined with symmetries of the external illumination may ease the solution of a radiative transfer problem significantly. For instance, a common assumption in spherical radiative transfer is to consider a model atmosphere which consists of homogeneous spherical shells. Furthermore, the solar source S 0 in Eq. (6) shows an axial symmetry with respect to the global zenith. It is due to both symmetries that the forward radiation field I does not depend on the global azimuth angle F [14,41]. Another well-known example with a high degree of symmetry is the plane-parallel radiative transfer problem assuming horizontally homogeneous layers with a uniform illumination at the top of the model atmosphere. Again, it is due to the symmetry of the model atmosphere and the symmetry of the radiation source that the forward radiation field shows a dependence only on height and on the direction of the radiation. Thus, the calculation of the forward radiation field is facilitated in both cases by virtue of the imposed symmetry relations. However, these simplifications cannot be adopted in a straight-forward manner in the calculation of the adjoint or pseudo-forward radiation field due to the fundamental difference between the solar source S0 and the response function R. In general, the response function has less symmetries and, in the case of the response function from Eq. (10), the full three-dimensional searchlight problem has to be solved, even for plane-parallel radiative transfer. Nevertheless, the symmetries of a radiative transfer problem can be used in a different

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

185

manner in order to simplify the formulation of the adjoint transfer problem. Here, we make use of the fact that the response function is not uniquely determined if the forward radiation field, as depicted above, is subject to symmetry properties. In other words, in cases in which the forward radiation field only depends on a reduced set of variables of the phase space, the original response function as well as its integration over the redundant variables describe the same radiative effect E. For example, the response function R in Eq. (10) and its integration Z 1 n dFRðr; Z; F; XÞ R ðr; Z; XÞ ¼ 2p 1 ¼ dðr  rTOA ÞdðZ  ZTOA ÞdðX  Xv Þ ð32Þ 2pr2 describe the same radiative effect E for the spherical shell atmosphere. Thus, in contrast to Eq. (8), the radiative effect E might be calculated via E ¼ hRn jIi Z Z dV ¼ V

dO Rn ðr; Z; XÞIðr; Z; XÞ:

ð33Þ

4p

The equivalent response function Rn represents, in comparison to the response function R from Eq. (10), an adjoint source which illuminates the top of the model atmosphere on a circle centered around the global zenith. This, in turn, introduces symmetries into the adjoint formulation analogous to those of the forward formulation. Therefore it suffices to perform adjoint and hence pseudo-forward calculations which obey the same simplified dependencies as the forward calculation. This corresponds to a reduction in dimensionality. Thus, Rn used as a pseudo-forward source, eases significantly the model calculations in the pseudo-forward formulation of radiative transfer, Eq. (24). Similar, in the case of plane-parallel radiative transfer, the associated pseudo-forward or search-light problem may be replaced by a corresponding radiative transfer problem with uniform illumination at the top of the model atmosphere. This leads to a concise formulation of the forward-adjoint perturbation theory in plane-parallel geometry [25,30]. Generally speaking, any symmetry of the forward radiative transfer problem can be mirrored to the adjoint formulation using the concept of an equivalent response function. This concept will be of use in the following section, in order to investigate a radiative transfer problem in simplified spherical geometry.

5. Application of the forward-adjoint perturbation theory to radiative transfer in a simplified spherical geometry 5.1. The radiative transfer model We study a simplified scenario in order to demonstrate the functionality of the forward-adjoint perturbation theory approach for a linearization of spherical radiative transfer. Here, the following approximations are made: (a) A spherical shell atmosphere is assumed, which means

ARTICLE IN PRESS 186

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

that the extinction and the scattering coefficient be and bs as well as the phase function P are constant within a spherical shell. (b) A radially symmetric illumination at the top of the atmosphere is assumed with m0 ¼ 1:0: Therefore, Eq. (6) can be expressed as S0 ðr; XÞ ¼ F 0 dðr  rTOA ÞdðjÞdðm þ 1Þ:

(34)

In this case the diffuse radiation field only depends on the radius r and the local zenith angle y: (See Fig. 2 for details). With these simplifications, the streaming term X  = in Eq. (3) can be written as X=¼m

q 1  m2 q þ ; qr r qm

(35)

which is the differential operator often used in astrophysical radiative transfer problems and neutron transport theory [12]. Due to the high degree of symmetry it is sufficient to solve the radiative transfer problem for the upper right quadrant only (Fig. 2). Despite these simplifications, the radiative transfer model still contains all important aspects of a spherical atmosphere. Furthermore, we assume a clear sky atmosphere, taking Rayleigh scattering and ozone absorption into account. Any refractional effects of radiation are neglected.

incoming radially symmetric radiation

r

θ

Ω

Fig. 2. The geometry of the simplified spherical radiative transfer model. The radially symmetric illumination at the top of the model atmosphere leads to a substantial simplification of the radiative transfer problem. Combined with a spherical shell atmosphere, the radiation field is independent from C; F and j: Thus, the radiation field only depends on the radius r and the local zenith angle y: Due to the high degree of symmetry, it is sufficient to calculate the radiation field in the upper right quadrant only.

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

187

5.2. The solution of the radiative transfer equation Using Eq. (35), the radiative transfer equation (3) can be written as   q 1  m2 q m þ Iðr; mÞ þ be ðrÞIðr; mÞ ¼ Sðr; mÞ: qr r qm

(36)

For the corresponding source function S from Eq. (4), Eq. (36) is a first order linear partial differential equation, which can be solved by the method of characteristics [42]. In Eq. (36), the streaming term describes the directional derivative along a characteristic line, which is given by a constant impact parameter [17,43,44] pffiffiffiffiffiffiffiffiffiffiffiffiffi (37) p ðr; mÞ ¼ r 1  m2 : Due to the fact that the impact parameter is the same for m; we have to distinguish between impact parameters for the up- and downward direction, which is indicated by an additional superscript þ and ; respectively. The slant optical depth along such a characteristic line p may be defined differentially as (38)

dt ¼ be ðrÞdðmrÞ;

where the product mr describes the geometrical pathlength. The new parameters provide a coordinate transformation from the ðr; mÞ-representation to a corresponding ðp ; tÞ-representation. In turn, Eq. (36) transforms into two ordinary differential equations, which may be integrated analytically over the slant optical depth t along a characteristic line p : With an additional splitting of the intensity field I in its direct and diffuse components I dir and I diff ; viz I ¼ I dir þ I diff ;

(39)

the resulting integral equations are given by  ðt2 t1 Þ I dir ðp  ; t1 Þ ¼ I dir ðp ; t2 Þe

(40)

and 



I diff ðp ; t1 Þ ¼ I diff ðp ; t2 Þe

ðt2 t1 Þ

Z

t2

þ

eðtt1 Þ Sdiff ðp ; tÞ dt;

(41)

t1

for two points at optical depth t1 and t2 on the corresponding characteristic line. In Eq. (40), p ¼ 0 indicates the characteristic line of the solar beam. It is I dir ðp  ; t ¼ 0Þ ¼ F 0 : The diffuse  scattering source function S diff at point ðp ; tÞ in Eq. (41) may be calculated by Z oðrÞ 1 Pðr; m0 ; mÞI diff ðr; m0 Þ dm0 Sðr; mÞ ¼ 2 1 Rr  be ðr0 Þ dr0 oðrÞ þ F 0 Pðr; m0 ; mÞe rTOA 4p Z 1 þ 2A Hðr; m0 ; mÞI diff ðr; m0 Þ dm0 1 Rr  be ðr0 Þ dr0 A þ F 0 m0 dðr  rBOA ÞYðmÞjmje rTOA ð42Þ p

ARTICLE IN PRESS 188

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

with oðrÞ ¼ bs ðrÞ=be ðrÞ: The first two terms in Eq. (42) describe the multiple scattering and the single scattering contribution to the source function. The last two terms in Eq. (42) describe the contribution of a reflecting surface to the source function. A subsequent transformation of Sðr; mÞ into the corresponding ðp ; tÞ-representation provides finally S diff ðp ; tÞ: For any numerical implementation of Eq. (41), the phase space in the ðp ; tÞ-representation has to be discretized in an appropriate manner. In this paper, we use an equidistant grid for p ; which is denoted by p i with i ¼ 0; . . . ; M: Furthermore, the model atmosphere is divided into K optically-homogeneous spherical shells. The intersection points of each characteristic line p i with the layer boundaries of the spherical shells provide a non-equidistant discretization tik with k ¼ 0; . . . ; K: It is tik1 ptotik : In particular, for p 0 ¼ 0 the parameter t00 corresponds to the vertical optical depth at the top of the atmosphere, whilst t0K corresponds to the vertical optical depth at the bottom of the atmosphere. Fig. 3 illustrates the chosen discretization. We use the renormalization scheme proposed by Grant and Hunt [45,46] in order to account for errors due to the discretization of the scattering phase function. Finally, this leads to 2ðM þ 1Þ coupled integral equations for the diffuse intensity. In order to perform the model calculations for the upper right quadrant only, it is necessary to take additional symmetry relations for characteristic lines in limb viewing geometry, i.e. p i XrBOA ; into account. At the tangent point, where the characteristic line enters or leaves the upper right quadrant (Fig. 3), the following relationship is valid  I diff ðp i ; tTP Þ ¼ I diff ðpi ; tTP Þ

(43)

with tTP the slant optical depth at the tangent point. For a further evaluation of the integrals in Eq. (41) two additional approximations are needed. First, the optically-homogeneous shells of the model atmosphere are split into optically-thin sublayers if the vertical optical thickness, thus for p i ¼ 0; exceeds a given threshold. Within the optically-thin sublayers the dependence of the intensity field on the optical depth t is approximated by its layer average. In principle, this introduces only a minor bias in the radiative transfer calculations depending on the threshold of the internal layer splitting. However, for p0

k=0

p1

.

.

.

p

M

UP

. . . k=K DOWN

Fig. 3. The upper right quadrant which is used for the calculation of the radiation field in Section 5 is discretized in ðM þ 1Þ ¼ 129 characteristic lines p i ; (i ¼ 0; . . . ; M) and K ¼ 100 homogeneous model layers. The radiation field is calculated along each characteristic line p i for upward (‘UP’) and downward (‘DOWN’) directed streams separately.

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

189

characteristic lines in limb viewing geometry, p i XrBOA ; this approximation becomes inaccurate in the lowest crossed atmospheric layer. Second, we replace the slant optical depth along the characteristic line p i in a spherical atmosphere by a Chapman function [16], viz tik  chapik :

(44)

Here, the Chapman functions are calculated at the grid points and any dependence on the pathlength is approximated by a linear interpolation between the boundary values. The 2ðM þ 1Þ coupled integral equations transform to a matrix equation with respect to the intensity field I diff with a numerical integration over m in Eq. (42) and due to the continuity of the intensity field at the layer boundaries [30]. This matrix equation can be solved numerically using standard computational techniques. In the following we employ the Gauss–Seidel iteration method [47], which was first proposed by Herman and Browning [48] to solve the radiative transfer equation in a plane-parallel atmosphere. It provides an efficient numerical scheme to solve the radiative transfer equation for clear sky and aerosol loaded atmospheres [30]. The Gauss–Seidel iteration method can be extended to spherical atmospheres in a straightforward way. 5.3. Verification of the simplified spherical radiative transfer model In order to verify the proposed radiative transfer model we compare calculated radiances at the top of the model atmosphere with those calculated by a Monte Carlo code. Details on the concept of Monte Carlo calculations and their application in radiative transfer simulations are given by e.g. Marchuk et al. [8]. The implementation of the simplified spherical geometry in a Monte Carlo code is straightforward and will not be discussed further. For all simulations we assume a model atmosphere of 100 km geometrical thickness, which has been divided into atmospheric layers of 1 km vertical extension. A total of 129 characteristic lines p i were used. Furthermore, we assume the satellite instrument to be located directly at the top of the model atmosphere. Fig. 4 shows a comparison between the proposed radiative transfer model and a Monte Carlo reference simulation for a surface albedo of A ¼ 0:0 and an Earth radius of RE ¼ 6378 km: Here, the radiance I TOA is normalized to the incoming solar flux F 0 : The upper graph in Fig. 4 shows the reflected radiance as a function of the zenith angle yv : A viewing geometry with yv ¼ 0:0o corresponds to a nadir observation of a satellite instrument, whereas yv 479:9o denote observations in limb viewing geometry. For yv ¼ 90:0o the satellite instrument looks horizontally into space and therefore would measure no radiation at all. The upper panel of Fig. 4 shows an increase of the radiance towards larger zenith angles, followed by a drop when the observation geometry reaches the limb. This spherical effect is in agreement with findings of Adams and Kattawar [10]. As yv increases, the length of the satellite’s line-of-sight also increases, which leads initially to a higher contribution of multiple-scattered photons to the radiance. However, extinction processes overweigh the multiple-scattering sources along the line-of-sight for larger zenith angles, leading to a decrease in the radiance. In limb viewing geometry the lineof-sight is so long that extinction processes dominate the radiative transfer. Additionally, the lineof-sight senses the upper part of the atmosphere, where less multiple scattering of photons takes place, as the zenith angle increases further. Both effects contribute to a drastic decrease in the radiance for limb viewing geometry.

ARTICLE IN PRESS 190

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

Fig. 4. Upper panel: Reflected radiance versus the viewing zenith angle yv : The radiance I TOA is normalized to the incoming radiative flux F 0 : The solid line shows the radiance calculated with the proposed spherical radiative transfer model (SPH) for a wavelength of 332 nm. The simulation is done for a clear-sky atmosphere including ozone absorption and Rayleigh scattering. The Earth radius is taken to be RE ¼ 6378 km: The surface albedo A is set to zero. In order to take the vertical inhomogeneity into account, a US-Standard atmosphere [49] is divided in 100 homogeneous layers of 1 km geometrical thickness. The dotted line shows the reflected radiance calculated with the spherical Monte Carlo reference code (MC). Lower panel: absolute difference DI TOA of the reflected radiance calculated with the proposed spherical radiative transfer model and the Monte Carlo reference code (solid line). For the Monte Carlo simulations pffiffiffiffiffiffiffiffiffia total number of N ¼ 107 photons were used. This results in an absolute statistical error which is proportional to N bin ; where N bin denotes the number of photons per discretized direction. The statistical error of the Monte Carlo calculations is depicted as a dashed line.

The lower panel in Fig. 4 denotes the absolute difference between the two radiance calculations. The comparison shows a slight underestimation of the reflected radiance calculated with the proposed spherical radiative transfer model. However, the calculated radiance stays within the bounds of the statistical error of the Monte Carlo simulation for most zenith angles. In limb viewing geometry, thus for yv 479:9o ; the underestimation of the reflected radiance becomes more serious and lies clearly outside the statistical error. This deviation is caused by the assumption of a linear behavior of the Chapman function, Eq. (44), within a model layer in the proposed radiative transfer model. This numerical approximation becomes inaccurate for the very large pathlengths in the lowermost atmospheric model layer of a limb viewing geometry. However, a finer discretization of this lowermost model layer may reduce the observed deviation significantly. A second comparison of the two radiative transfer models has been carried out for a Lambertian surface albedo of A ¼ 1 (Fig. 5). The upper panel of Fig. 5 once again shows the reflected radiances with respect to the zenith angle. Due to the reflecting surface the intensity is

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

191

Fig. 5. Same as Fig. 4 but for a surface albedo of A ¼ 1:

higher compared to the case depicted in Fig. 4. Moreover, the radiance at the top of the atmosphere decreases continuously with increasing zenith angle. The radiation field is dominated by surface effects and an increase of the line-of-sight leads directly to a decrease in the reflected radiance due to an enhanced extinction. The lower panel of Fig. 5 again shows the differences in the reflected radiance between the proposed radiative transfer model and a Monte Carlo reference simulation. Both models agree very well for almost all viewing geometries. An underestimation of the reflected radiance calculated with the proposed spherical radiative transfer model can once again be seen for yv 479:9o : The reasons are the same as discussed before. A third comparison has been performed for a planetary radius of RE ¼ 63:78 km with a surface albedo of A ¼ 0:0: The upper panel of Fig. 6 shows again the reflected radiances with respect to the zenith angle and the lower panel of Fig. 6 demonstrates the accuracy of the proposed radiative transfer model with respect to the Monte Carlo reference simulation. The proposed model slightly overestimates the result of the Monte Carlo simulation for small zenith angles. Nevertheless, the general agreement is good. Limb viewing geometry already occurs for zenith angles yv 422:9 in the case of this very small planetary radius. In order to investigate the effect of the curvature of the Earth’s atmosphere on the radiation field we compared the reflected radiance calculated in spherical geometry I sph with the radiance calculated in plane-parallel geometry I pp : Fig. 7 shows the ratio of the two quantities. The planetary radius changes between 6:378 km and 6:378  1013 km for the spherical simulations. It is expected that the spherical radiative transfer model behaves like a plane-parallel one for the very large radius of RE ¼ 6:378  1013 km: In this case the ratio I sph =I pp almost reaches unity for all zenith angles. Also for the Earth radius RE ¼ 6378:0 km the reflected radiance in spherical

ARTICLE IN PRESS 192

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

Fig. 6. Same as Fig. 4 but for a planetary radius of RE ¼ 63:78 km:

Fig. 7. Ratios of the reflected radiance calculated for the simplified spherical geometry, I sph ; and the reflected radiance calculated in plane-parallel approximation, I pp ; versus the viewing zenith angle yv : Same atmospheric conditions as described in the caption of Fig. 4. The planetary radii range from RE ¼ 6:378 km to RE ¼ 6:378  1013 km:

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

193

geometry differs insignificantly from a corresponding plane-parallel calculation for zenith angles yv o75 : However, for zenith angles larger than 75 ; spherical effects cause a clear drop in the ratio I sph =I pp : The spherical shape of the atmosphere becomes more and more important for even smaller planetary radii (RE ¼ 637:8 km; RE ¼ 63:78 km and RE ¼ 6:738 km). In these cases spherical effects occur even for small and intermediate zenith angles and the reflected radiances show distinctive peaks. In nadir viewing geometry the ratio I sph =I pp decreases with smaller planetary radius. This spherical effect is caused by the increase of the photon escape probability at the top of the model atmosphere for a stronger curvature of the atmosphere. With the verification conducted in this section, we conclude that the proposed spherical radiative transfer model is well-suited for calculating the radiation field in a simplified spherical atmosphere with sufficient accuracy. 5.4. Solution of the pseudo-forward transfer equation In order to calculate the derivatives of the reflected radiance using Eq. (31), one has to solve the radiative transfer problem in its forward and pseudo-forward formulation, Eq. (12) and Eq. (24), respectively, for the simplified spherical geometry. In order to solve the pseudo-forward transfer equation (24), we consider a nadir-viewing satellite instrument. For this particular viewing geometry, as is discussed in the previous section, spherical effects become important for small planetary radii. Due to the symmetries of the forward radiative transfer problem, the response function in Eq. (10) has an equivalent response function 1 dðr  rTOA Þdðm  1Þ (45) 8p2 r2 after integration over the cosine of the global zenith angle Z; the global azimuth angle F and the local azimuth angle j (see Section 4). The searchlight problem simplifies therefore to a radiative transfer problem with a uniform pseudo-forward illumination at the top of the model atmosphere. Hence, the forward and the pseudo-forward problem are identical except for a different strength of the external sources S 0 and S x in Eqs. (12) and (24), respectively. Therefore the pseudo-forward radiance field is given by Rn ðr; XÞ ¼

1 Iðr; XÞ: (46) 8p2 r2 Thus, we obtain the solution of the corresponding pseudo-forward problem x; multiplied with the factor 1=ð8p2 r2 Þ; simultaneously with the solution of the forward radiative transfer problem. This simplification is only possible for a nadir viewing scenario, which makes it useful for demonstration purposes. xðr; XÞ ¼

5.5. Linearization of the simplified spherical radiative transfer model For a given forward and pseudo-forward radiance field, Eq. (31) provides an analytical approach for calculating the derivative of the reflected radiance with respect to the trace gas number density x in a model layer k. Due to the symmetries of the forward and pseudo-forward problem for the simplified geometry and for a nadir-viewing satellite instrument,

ARTICLE IN PRESS 194

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

Eq. (31) reduces to sk qI TOA ðx0 Þ ¼ abs qx 2

Z

rk

Z

1

dr dm I 0 ðr; mÞI 0 ðr; mÞ; rk1

(47)

1

where we have made use of the fact that the radiative effect E can be identified with the radiance at the top of the model atmosphere I TOA : The radiance field I 0 ðr; mÞ is the solution of the forward radiative transfer problem for the unperturbed atmospheric parameter x0 in Eq. (47). Furthermore, the absorption cross-section of the considered trace gas species skabs is assumed to be constant within a model layer k. The perturbation integrals in Eq. (47) can be evaluated numerically with low computational costs. In order to verify the proposed approach, we choose ozone as the relevant trace gas species to be investigated. For reference purposes we use the finite difference approach I TOA ðx0 þ DxÞ  I TOA ðx0 Þ qI TOA ! for Dx ! 0; (48) Dx qx x0

Fig. 8. Left panel: derivative of the reflected radiance I TOA normalized to the incoming solar flux with respect to the ozone density at different altitude levels. The solid line denotes the derivative calculated with the spherical perturbation theory approach (SPT) from Eq. (47). The dashed line denotes the derivative calculated with the finite difference scheme (FD), Eq. (48). The computations are performed at a wavelength of 332 nm for clear sky conditions (see Fig. 4). The radius of the planet is taken to be RE ¼ 6378 km: A relative perturbation in the ozone number density of Dx ¼ 102 x throughout the atmosphere is sufficient for a reliable determination of qI TOA =qx with the finite difference scheme. Right panel: relative difference d between computations using the spherical perturbation theory and the finite difference scheme.

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

195

which converges to the derivative for small perturbations in the ozone number density Dx: The perturbed and the unperturbed radiance at the top of the model atmosphere, I TOA ðx0 þ DxÞ and I TOA ðx0 Þ; are calculated with the forward mode of the simplified radiative transfer model, which is also used for the evaluation of the perturbation integrals in Eq. (47). Fig. 8 shows the derivative qI TOA =qx calculated with the spherical perturbation theory approach and the finite difference approach as well as the relative deviation between the two methods for a wavelength of 332 nm and a planetary radius of RE ¼ 6378 km: The calculation of the derivatives is performed for nadir viewing geometry with yv ¼ 0:0o and for a surface albedo of A ¼ 0: The sensitivity of I TOA with respect to changes in the ozone number density shows a weak dependence on altitude with a maximum at around 50 km (panel on the left-hand side of Fig. 8). Below 5 km altitude the sensitivity decreases to zero, which is due to the vanishing surface reflection. The relative difference between the perturbation theory approach and the finite difference approach is depicted on the panel on the right-hand side of Fig. 8. The overall deviation is less than 0.2%. Both approaches agree very well in the altitude range from 30 to 100 km. Slight deviations occur below an altitude of 30 km, which are caused by the layer-averaging of the multiple-scattered radiance field (compare Section 5.2). This approximation becomes inaccurate for the tangent layer of a characteristic line in limb viewing geometry. A second comparison was carried out for a planetary radius of RE ¼ 63:78 km: For such a radius, significant spherical effects occur in nadir viewing geometry (compare Fig. 7). The panel on the left-hand side of Fig. 9 once again shows the derivatives qI TOA =qx calculated with the spherical perturbation theory approach and the finite difference scheme. The sensitivity throughout the atmosphere is slightly smaller than in the previous case. This may be due to the

Fig. 9. Same as Fig. 8 but for a planetary radius of RE ¼ 63:78 km:

ARTICLE IN PRESS 196

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

enhanced photon escape probability for smaller planetary radii, leading to a decreased interaction of photons with the participating medium. The vertical shape of the sensitivity, however, stays approximately the same. The relative difference of the perturbation theory approach to the finite difference approach is depicted in the panel on the right-hand side of Fig. 9. Here, the overall deviation is smaller than 1.0%, with the exception of the lowest atmospheric layer where the relative deviation reaches an absolute value of 3.1%. However, as the sensitivity in this atmospheric region is very small for ozone, this deviation is of minor importance. The deviations, which occur in the lower atmosphere, again are caused by the layer-averaging of the multiplescattering radiance field. These deviations are slightly larger than in the previous case, as more characteristic lines are in limb viewing geometry due to the very strong bending of the atmosphere associated with the small radius. Despite the observed deviations in the lower atmosphere, we can conclude that the spherical perturbation theory approach is suited for an accurate calculation of the derivatives qI TOA =qx:

6. Conclusions We have presented a general linearization approach for radiative transfer in spherical geometry using the forward-adjoint perturbation theory. For this purpose the radiation transport in a spherical atmosphere is described via an operator formalism. This enables the derivation of the corresponding linearization principles in a universally valid notation. The expressions for the derivatives presented in this paper can be applied to linearize any spherical radiative transfer code that calculates the internal forward and adjoint radiation field. The essential advantage of the spherical perturbation theory with respect to the finite difference approach is its computational efficiency. The finite difference approach requires a new computation of the intensity field for each perturbation of an atmospheric parameter. The perturbation theory approach, however, enables one to calculate the derivatives with respect to this parameter at any altitude based on the forward and adjoint intensity field only. Thus, the spherical perturbation theory approach provides a highly efficient linearization of spherical radiative transfer. In order to demonstrate the capability of this method we developed a linearized spherical radiative transfer model for a simplified spherical atmosphere. The assumption of homogeneous spherical shells and a radially symmetric illumination of the sphere leads to radiative transfer calculations, which only depend on the radius and the local zenith angle. The radiative transfer model presented here provides the radiation field as well as the derivatives of the reflected radiance with respect to atmospheric trace gas profiles. The radiances calculated with the proposed model are verified with a Monte Carlo simulation in spherical geometry for different atmospheric conditions and planetary radii. At a wavelength of l ¼ 332 nm the deviation of the spherical radiative transfer model from the Monte Carlo reference model is less than approximately 0.5% for zenith angles smaller than 80 : Deviations occur for zenith angles larger than 80 ; which are due to the linear interpolation of the slant optical thickness, expressed by a Chapman function, in a model layer in the proposed radiative transfer model. Furthermore, the accuracy of the derivatives of the reflected radiance with respect to the ozone density was investigated for a nadir viewing scenario. Differences between the derivatives calculated with the perturbation theory approach and numerical calculations using a finite difference scheme are less than 0.2%

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

197

throughout the atmosphere for a planetary radius of RE ¼ 6378 km; and less than 1.0% for a radius of RE ¼ 63:78 km: As the basic principles used for the linearization of the radiative transfer model do not only hold for a simplified geometry but are formulated quite generally, we conclude that the presented forward-adjoint perturbation theory approach can be used for an efficient linearization of fully spherical radiative transfer, which is currently a subject of investigation.

Acknowledgments The authors acknowledge Marco Spaans and Francesco Spada for useful discussions on the Monte Carlo code. We would like to thank Ahilleas Maurellis for comments and suggestions on an earlier version of this paper. This research was supported by SRON under project number 6430-SCIARALI (GO-2).

Appendix A. The adjoint transport operator In this section it is shown that Ly from Eq. (16) is indeed the adjoint operator which fulfills Eq. (14). Here, we follow the outline presented by Bell and Glasstone [23]. First we redefine the boundary conditions for the forward and the adjoint intensity, Eqs. (7) and (15), as IðrTOA ; XÞ ¼ 0

for  1pn  Xo0;

IðrBOA ; XÞ ¼ 0

for  1pn  Xo0

I y ðrTOA ; XÞ ¼ 0

for 0pn  Xo1;

ðA:1Þ

and y

I ðrBOA ; XÞ ¼ 0

for 0pn  Xo1;

ðA:2Þ

respectively. Here, n is a unit vector in the direction of the outward normal at a position r on the boundary. Any radiation with n  X40 crosses the surface in an outward direction whereas radiation for which n  Xo0 crosses the surface in an inward direction. Furthermore, Eq. (14) can be rewritten as Z Z Z Z y dV dO I LI ¼ dV dO ILy I y : (A.3) V

V

4p

4p

It can be shown in a straightforward manner that the difference between the right and the lefthand side of Eq. (A.3) is given by the difference of the terms, which contain the streaming terms of I y LI and ILy I y ; viz Z Z   D¼ dV dO I y ðX  =IÞ þ IðX  =I y Þ : (A.4) V

4p

In the following it is demonstrated that this difference also vanishes, thus D ¼ 0:

ARTICLE IN PRESS 198

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

Since the nabla operator = does not operate on directions, it is allowed to write X  =I ¼ =  XI

(A.5)

X  =I y ¼ =  XI y :

(A.6)

and

Then the two terms in Eq. (A.4) may be combined as Z Z D¼ dV dO =  XI y I: V

(A.7)

4p

The first integral can be converted to a surface integral by the use of the divergence theorem, which yields Z Z D¼ dA dO n  XI y ðr; XÞIðr; XÞ; (A.8) A

4p

where the surface integration is performed over the bounding surface A, on which the boundary conditions are imposed. On this surface, the product I y I is always zero. According to the boundary conditions, Eqs. (A.1) and (A.2), I y is zero for n  XX0 and I is zero for n  Xo0: Consequently, D ¼ 0 and hence Eq. (14) is satisfied.

References [1] Bovensmann H, Burrows JP, Buchwitz M, Frerick J, Noel S, Rozanov VV, Chance KV, Goede APH. Sciamachy: mission objectives and measurement modes. J Atmos Sci 1999;56:127–50. [2] van der A RJ. Improved ozone profile retrieval from combined nadir/limb observations of sciamachy. J Geophys Res 2001;106:14583–94. [3] Payne WF, Llewellyn EJ, Matsushita JS. Osiris: an imaging spectrograph for odin. 46th International Astronautical Congress, Oslo, Norway; 1995. [4] Phillips P. A technique for the numerical solution of certain integral equations of the first kind. J Assoc Comput Mach 1962;9:84–97. [5] Tikhonov A. On the solution of incorrectly stated problems and a method of regularization. Dokl Akad Nauk SSSR 1963;151:501–4. [6] Rodgers C. Retrieval of atmospheric temperature and composition from remote measurements of thermal radiation. Rev Geophys 1976;14:609–24. [7] Hasekamp O, Landgraf J. Ozone profile retrieval from backscattered ultraviolet radiances: the inverse problem solved by regularization. J Geophys Res 2001;106:8077–88. [8] Marchuk GI, Mikhailov GA, Nazaraliev MA, Darbinjan RA, Kargin BA, Elepov BS. The Monte Carlo methods in atmospheric optics. Berlin: Springer; 1980. [9] Collins DG, Blaettner WG, Wells MB, Horak HG. Backward monte carlo calculations of the polarization characteristics of the radiation emerging from spherical shell atmospheres. Appl Opt 1972;11:2684–96. [10] Adams CN, Kattawar GW. Radiative transfer in spherical shel atmospheres, i. rayleigh scattering. Icarus 1978;35:139–51. [11] Oikarinen L, Sihvola E, Kyrola E. Multiple scattering in limb viewing geometry. J Geophys Res 1999;104:31261–74. [12] Sen KK, Wilson SJ. Radiative transfer in curved media. Singapore: World Scientific Publishing; 1990. [13] Sobolev VV, Minin IN. Light scattering in the spherical atmosphere. ISZ (artificial earth satellites) 14, 1962. [14] Sobolev V. Light scattering in planetary atmospheres. Oxford: Pergamon Press; 1975.

ARTICLE IN PRESS HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

199

[15] Smokty OI. Multiple light scattering in the spherical planetary atmosphere. Pure Appl Geophys 1969;72:214–26. [16] Dahlback A, Stamnes K. A new spherical model for computing the radiation field available for photolysis and heating at twilight. Planet Space Sci 1991;39:671–83. [17] Lary DJ, Balluch M. Solar heating rates: the importance of spherical geometry. J Atmos Sci 1993;50:3983–93. [18] Balluch M. A new numerical model to comput photolysis rates and solar heating with anisotropic scattering in spherical geometry. Ann Geophys 1996;14:80–97. [19] Anderson DE. The troposphere–stratosphere radiation field at twilight: a spherical model. Planet Space Sci 1983;31:1517–23. [20] Herman BM, Ben-David A, Thome KJ. Numerical technique for solving the radiative transfer equation for a spherical shell atmosphere. Appl Opt 1994;33:1760–70. [21] Rozanov A, Rozanov V, Burrows JP. A numerical radiative transfer model for a spherical planetary atmosphere: combined differential-integral approach involving the picard iterative approximation. JQSRT 2001;69:491–512. [22] Kaiser JW, Burrows JP. Fast weighting functions for retrievals from limb scattering measurements. JQSRT 2003;77:273–83. [23] Bell GI, Glasstone S. Nuclear reactor theory. New York: Van Nostrand Reinholt; 1970. [24] Marchuk G. Equation for the value of information from weather satellites and formulation of inverse problems. Cosmic Res 1964;2:394–409. [25] Box M, Gerstl S, Simmer C. Application of the adjoint formulation to the calculation of atmospheric radiative effects. Beitr Phys Atmos 1988;61:303–11. [26] Ustinov EA. Inverse problem of the photometry of the solar radiation reflected by an optically thick planetary atmosphere, mathematical framework and weighting functions of linearized inverse problem. Cosmic Res 1991;29:519–31. [27] Ustinov EA. Inverse problem of the photometry of the solar radiation reflected by an opticaly thick planetary atmosphere 2. Numerical aspects and requirements on the observation geometry. Cosmic Res 1991;29:785–800. [28] Ustinov EA. Inverse problem of the photometry of the solar radiation reflected by an optically thick planetary atmosphere 3. Remote sensing of minor gaseous constituents and an atmospheric aerosol. Cosmic Res 1992;30:170–81. [29] Rozanov VV, Kurosu T, Burrows JP. Retrieval of atmospheric constituents in the uv-visible: a new quasianalytical approach for the calculation of weighting functions. JQSRT 1998;60:277–99. [30] Landgraf J, Hasekamp OP, Trautmann T, Box MA. A linearized radiative transfer model for ozone profile retrieval using the analytical forward-adjoint perturbation theory approach. J Geophys Res 2001;106:27291–306. [31] Hasekamp OP, Landgraf J, van Oss R. The need of polarization modeling for ozone profile retrieval from backscattered sunlight. J Geophys Res 2002;107:4692. [32] Chandrasekhar S. Radiative transfer. New York: Dover Publications, Inc.; 1960. [33] Mishchenko MI, Lacis AA, Travis LD. Errors induced by the neglect of polarization in radiance calculations for rayleigh-scattering atmospheres. JQSRT 1994;51:491–510. [34] Landgraf J, Hasekamp OP, Trautmann T. Linearization of radiative transfer with respect to surface properties. JQSRT 2002;72:327–39. [35] Lewins J. Importance. The adjoint function. The physical basis of variational and perturbation theory in transport and diffusion problems. New York: Pergamon Press; 1965. [36] Case KM, Zweifel PF. Linear transport theory. Reading: Addison-Wesley Publishing Co; 1967. [37] Case KM. Transfer problems and the reciprocity principle. Rev Mod Phys 1957;29:651–63. [38] Chandrasekhar S. On the diffuse reflection of a pencil of radiation by a plane-parallel atmosphere. Proc Natl Acad Sci 1958;44:933–40. [39] Box M, Gerstl S, Simmer C. Computation of atmospheric radiative effects via perturbation theory. Beitr Phys Atmos 1989;62:193–9. [40] Walter HH, Landgraf J, Hasekamp OP. Linearization of a pseudo-spherical vector radiative transfer model. JQSRT 2004;85:251–83. [41] Lenoble J, Sekera Z. Equation of radiative transfer in a planetary spherical atmosphere. Proc Natl Acad Sci USA 1961;47:372–89. [42] Bleecker D, Csordas G. Basic partial differential equations. International Press Publications; 1996.

ARTICLE IN PRESS 200

HH. Walter, J. Landgraf / Journal of Quantitative Spectroscopy & Radiative Transfer 95 (2005) 175–200

[43] Hummer DG, Rybicki GB. Radiative transfer in spherically symmetric systems, the conservative grey case. Mon Not R Astron Soc 1971;152:1–19. [44] Yorke HW. Numerical solution of the equation of radiation transfer in spherical geometry. Astron Astrophys 1980;86:286–94. [45] Grant IR, Hunt GE. Discrete space theory of radiative transfer: I, fundamentals. Proc R Soc London A 1969;313:183–97. [46] Wiscombe WJ. On initialization, error and flux conservation in the doubling method. JQSRT 1976;16:637–8. [47] Strang G. Linear algebra and its applications. San Diego, CA: Harcourt Brace Jovanovich; 1986. [48] Herman B, Browning S. A numerical solution to the equation of radiative transfer. J Atmos Sci 1965;22:559–66. [49] National Oceanic and Atmospheric Administration (NOAA), U.S. Standard Atmosphere, Rep. NOAA-S/ T76-1562, Washington, DC, US Government Printing Office, 1976.