Zdunkowski W., Trautmann T., Bott A. Radiation in the atmosphere: A course in Theoretical Meteorology
Подождите немного. Документ загружается.
4.6 The Monte Carlo method 125
A new random number R
t,1
will determine the path length s
1
between P
1
and
the new point P
2
. Substituting this path length together with (4.137) into (4.127)
yields the new coordinates (x
2
, y
2
, z
2
) of the point P
2
. At this point new local
scattering angles are chosen as random numbers and the computational process is
repeated until the test photon either is absorbed within the medium or at the ground
or it leaves the atmosphere into space, (event E in Figure 4.10). Since reflection at
the ground is assumed to be isotropic, two uniformly distributed random numbers
from the intervals [0,π/2) and [0, 2π) can be used to find the new direction of the
test photon after ground reflection.
The above discussion was based on the assumption that the medium is homoge-
neous. In case of a vertically and horizontally inhomogeneous medium the transmis-
sion calculation must fully account for the (x, y, z)-dependence when determining
the optical depth, i.e. (4.128) has to be generalized to
T (s) = exp
−
s
k
ext
(s) ds
(4.138)
In addition to the already defined levels z
j
for the vertically inhomogeneous
medium, we have to introduce similar grids for discretizing the (x, y)-space. There-
fore, for a general three-dimensional medium the space is discretized into small
volume elements V . The photon paths are then traced through these individ-
ual volume elements. Radiative fluxes or actinic fluxes may be determined at the
midpoints of the six faces of each volume element by weighting them with the
corresponding projection factor as the photons intersect a particular reference area.
In case of actinic fluxes this projection factor is always 1 since the weight of
each test photon is independent of its flight direction. For radiative flux densities
through area elements x = con st the weighting factor is equal to the cosine of
the angle α subtended by the outward normal n
x
of the area element and the unit
vector Ω specifying the flight direction of the photon, see Figure 4.13. The energy
is counted positive if the photon travels into the interior of V , otherwise it is
negative.
4.6.2 Treatment of absorption
The initial energy of the model photons is reduced by absorption due to gases and
atmospheric particles. Let N represent the total number of model photons, e.g.
some 10 000, entering the top of the atmosphere. On entry the initial energy carried
by such a photon is given by E
0
= µ
0
S
0
/N where S
0
refers to the solar constant
within a small spectral interval. Let s
i, j
specify the total photon path between the
starting point P
0
and an arbitrary intersection point D
i, j
as defined in Figure 4.10.
126 Quasi-exact solution methods for the RTE
x
y
z
Ω
n
x
α
Fig. 4.13 Definition of the angle α between the photon’s flight direction Ω and
the outward normal n
x
of the area x = const of the volume element V .
The length of this path is then given by
s
i, j
=
i−1
k=0
s
k
+ s
i, j
(4.139)
where s
i, j
is the distance between point P
i
and D
i, j
.Ifω
0,i
is the single scattering
albedo of the particulate material located at point P
i
, then the energy of the model
photon after traversing the path s
i
, may be expressed as
E(s
i, j
) = E
0
[1 − A(s
i, j
)]
i
!
k=1
ω
0,k
(4.140)
where A(s
i, j
) is the so-called absorption function of the atmospheric absorber gas.
This function will be discussed in detail in a later chapter. In the special case of
reflection at the ground the value of the corresponding ω
0,k
has to be replaced
by the ground albedo A
g
. From (4.140) it is seen that the effect of gaseous and
particulate absorption can be determined after the photon trajectories have been
computed in a purely scattering atmosphere. Some authors treat the absorption
process differently, namely by reducing the photon’s energy at each interaction
point during the simulation of each photon trajectory.
Finally, we have to decide what kind of interaction takes place at a particular inter-
action point P
i
. In general, three different processes are possible, that is molecular
scattering, extinction by aerosol particles or extinction by cloud droplets. Gaseous
absorption is already treated in terms of the absorption function A(s
i, j
) in (4.140).
In order to determine the kind of interaction at P
i
additional random numbers are
drawn. Usually it is assumed that only one of the possible interactions takes place.
However, some authors also allow for combined extinction processes of aerosols
and water drops. In these cases average values of the extinction ccefficients and
4.6 The Monte Carlo method 127
phase functions of the different particles have to be calculated. Monte Carlo algo-
rithms designed in the above manner have been employed, for example, by O’Hirok
and Gautier (1998a,b) and by Trautmann et al. (1999), to name a few. Radiative
characteristics of finite cloud fields have also been calculated using MCM, see e.g.
Welch and Zdunkowski (1981). For more details the interested reader is referred to
the original literature.
In the previous subsections we described how the random paths of test photons
can be determined. In order to obtain a statistically sound result a huge number N of
such photon trajectories has to be simulated. All we need to do for the computation
of the radiance field and the flux density is to count the number and energy of the test
photons passing all intersection points D
j
=
i
D
i, j
at height z
j
. For simplicity
we assume a horizontally homogeneous atmosphere so that all photon paths can
be related to a single vertical atmospheric reference column having unit area. This
treatment can be justified by the fact that for uniform external illumination in each
such reference column, on average, the same number of photons with a certain
energy and flight direction can be found.
Let us assume that for D
j
a total of J (z
j
) test photons are registered. Let us
further introduce a total of 2s solid angle elements
i
which cover the 4π unit
sphere, i.e.
lower 2π hemisphere:
i
, i = 1,...,s
upper 2π hemisphere:
i
, i = s + 1,...,2s
(4.141)
The number of test photons confined in
i
is designated by J
i
(z
j
) so that
J (z
j
) =
2s
i=1
J
i
(z
j
) (4.142)
Each of the J
i
(z
j
) photons may carry a particular energy E
l
(
i
, z
j
) which is the
energy of model photon number l within
i
penetrating the reference surface z
j
.
Adding the energy contributions of all photons yields the radiance at level z
j
within
the solid angle element
i
I
i
(z
j
) =
1
i
tF
J
i
(z
j
)
l=1
E
l
(
i
, z
j
)
(4.143)
Here t and F are the time interval and the reference area referring to photons
entering the top of the atmosphere at point P
0
. From (4.143) we may compute the
128 Quasi-exact solution methods for the RTE
upward and downward flux densities as
E
+
(z
j
) =
2s
i=s+1
I
i
(z
j
)
i
, E
−
(z
j
) =
s
i=1
I
i
(z
j
)
i
(4.144)
In a three-dimensional atmosphere similar expressions can be defined to determine
the flux densities crossing area elements x = const and y = const. We will not
discuss this topic any further.
The Monte Carlo method described in the previous subsections is not limited to
a horizontally homogeneous situation and for media infinitely extending in the lat-
eral directions. Obeying the general principles introduced above, one can construct
a Monte Carlo model for a three-dimensional situation. This can be achieved by
introducing in analogy to the area z
j
= con st vertically oriented areas x = const
and y = const with corresponding arrays to register the intersections of the pho-
tons through these area elements. Moreover, the optical properties for the extinction
coefficient, the single scattering albedo and the details of the scattering phase func-
tion must be supplied to represent the physical characteristics of the finite medium.
In addition, one has to realize that solar photons not only illuminate the upper
boundary of the model domain but also enter it through the lateral sunlit surfaces.
A comprehensive description of the MCM may be found in Davis et al. (1979).
A full discussion of the rotation matrices determining the flight direction of a model
photon after a scattering event is given in Zdunkowski and Korb (1985). For various
applications of the MCM to two- and three-dimensional spatially inhomogeneous
media we refer, for example, to the work of Barker and Davies (1992), Cahalan
et al. (1994), Los et al. (1997), and more recently to Trautmann et al. (1999).
4.7 Appendix
4.7.1 The reflection matrix at the ground
We will now derive the reflection matrix r
m
g
for m = 0 under the assumption of
isotropic ground reflection for both the incident diffuse light and the direct beam.
In addition, we include possible contributions due to the thermal radiation of the
ground.
The upward directed flux density at the ground (τ = τ
N
)isgivenby
E
+,z
(τ
N
) = A
g
E
−,z
(τ
N
) + µ
0
S
0
exp
−
τ
N
µ
0
+ (1 − A
g
)π B
g
(4.145)
where A
g
is the albedo of the isotropically reflecting ground, E
±,z
(τ
N
) are the
upward and downward traveling flux densities at ground level, µ
0
S
0
exp(−τ
N
/µ
0
)
4.7 Appendix 129
the direct solar flux density reaching the ground, and B
g
the black body emission
of the ground with temperature T
g
.
Owing to the isotropy of the reflected radiation the flux densities in (4.145) can
be replaced by corresponding isotropic radiances, i.e.
2π
1
0
µI
m=0
+
(τ
N
)dµ = A
g
2π
1
0
µI
m=0
−
(τ
N
,µ)dµ + A
g
µ
0
S
0
exp
−
τ
N
µ
0
δ
0m
+(1 − A
g
)π B
g
δ
0m
(4.146)
The integral on the left-hand side of this equation may be evaluated yielding
2π
1
0
µI
m=0
+
(τ
N
)dµ = 2π I
m=0
+
(τ
N
)
µ
2
2
µ=1
µ=0
= π I
m=0
+
(τ
N
) (4.147)
Expressing in (4.146) the integral containing the downwelling radiation by means
of the Gaussian quadrature gives together with (4.147)
I
m=0
+
(τ
N
) = 2A
g
s
i=1
w
i
µ
i
I
m=0
−
(τ
N
,µ
i
) +
A
g
π
µ
0
S
0
exp
−
τ
N
µ
0
δ
0m
+(1 − A
g
)B
g
δ
0m
(4.148)
Since the reflected radiation is isotropic it follows that I
m=0
+
(τ
N
) has the same value
for each direction µ
i
(i = 1,...,s). Therefore, the upwelling radiation can also be
written in vector–matrix form as
I
m=0
+
(τ
N
)
.
.
.
I
m=0
+
(τ
N
)
= 2A
g
µ
1
w
1
... µ
s
w
s
.
.
.
.
.
.
.
.
.
µ
1
w
1
... µ
s
w
s
I
m=0
−
(τ
N
,µ
1
)
.
.
.
I
m=0
−
(τ
N
,µ
s
)
+
A
g
π
µ
0
S
0
exp
−
τ
N
µ
0
δ
0m
1
.
.
.
1
+ (1 − A
g
)B
g
δ
0m
1
.
.
.
1
(4.149)
Substitution of the source vectors J
m=0
+,1
(N + 1, N ), J
m=0
+,2
(N + 1, N ) from (4.17)
leads to the final expression for the reflected radiance
I
m=0
+
(τ
N
) = 2A
g
µ
1
w
1
... µ
s
w
s
.
.
.
.
.
.
.
.
.
µ
1
w
1
... µ
s
w
s
I
m=0
−
(τ
N
)
+ J
m=0
+,1
(N + 1, N ) +J
m=0
+,2
(N + 1, N ) (4.150)
130 Quasi-exact solution methods for the RTE
At τ = τ
N
the upwelling radiation may be written as
I
m=0
+
(τ
N
) = r
m=0
g
I
m=0
−
(τ
N
) + J
m=0
+,1
(N + 1, N ) +J
m=0
+,2
(N + 1, N ) (4.151)
Comparison of this expression with (4.150) gives the explicit form of the reflection
matrix
r
m=0
g
= 2 A
g
µ
1
w
1
... µ
s
w
s
.
.
.
.
.
.
.
.
.
µ
1
w
1
... µ
s
w
s
(4.152)
as stated in (4.18).
It should be mentioned that for anisotropic ground reflection a similar form of
the reflection matrix has to be derived from the reflection function of the ground for
each azimuthal expansion term m > 0. In this context it is interesting to note that
for an isotropically reflecting ground the boundary conditions for the azimuthal
expansion terms of the radiance field I
m≥1
+
(τ
N
) are identical with the so-called
vacuum boundary conditions. Vacuum boundary conditions mean that there exists
no diffuse radiation that would enter the medium at its upper or lower boundary.
The treatment of complex reflection laws at the ground, that is the spec-
ification of the so-called bidirectional reflection function for diffuse radiation
ρ
d
(µ, ϕ, −µ
,ϕ
), is rather difficult. Some information for treating anisotropic
ground reflection within the radiative transfer problem can be found in Stamnes
et al. (1988). These authors assume that ρ
d
is a function only of the angle between
the directions of incidence and reflection, i.e. there is no preferred direction for
the scattering of radiation at the lower boundary. In this case, and in analogy to
the phase function, the bidirectional reflection function can be expanded into its
Fourier components
ρ
d
(µ, ϕ, −µ
,ϕ
) =
2n−1
m=0
(2 − δ
0m
)ρ
m
d
(µ, −µ
) cos m(ϕ − ϕ
) (4.153)
The contribution of the ground to the upwelling radiance field due to diffuse
anisotropic reflection of the downwelling diffuse and direct radiation can then be
expressed by
I
m
(µ, τ
N
) = (1 + δ
0m
)
1
0
ρ
m
d
(µ, −µ
)µ
I
m
(τ
N
, −µ
) dµ
+µ
0
S
0
exp
−
τ
N
µ
0
ρ
m
d
(µ, −µ
0
) (4.154)
4.8 Problems 131
4.8 Problems
4.1: To become familiar with the terminology of the addition theorems and the
interaction principle check the validity of equations (4.4) and (4.5).
4.2: Check in detail equations (4.7) and (4.8) by carrying out the comparison.
Note well that this problem is not designed to keep you busy (as it may
appear) but to familiarize yourself with the general procedure of the matrix
operator method.
4.3: Carry out in detail the comparison to verify (4.13).
4.4: Consider an inhomogeneous atmosphere consisting of three homogeneous
sublayers. Each sublayer has different optical properties. Follow the out-
lined procedure that summarizes the main steps of the MOM to formally
obtain the vectors I
m
+
(τ = 0) and I
m
±
(τ
3
) and the upward and downward radi-
ance vectors at the sublayer boundaries. Make sure you have the required
information before you continue with the next step.
4.5: Since Gaussian quadrature (Gauss–Legendre quadrature) is essential in our
work, you should get an impression of the accuracy of this method. For
this purpose evaluate the following integral I =
1.5
0.2
exp(−x
2
)dx by using
a three-term formula. The weights w
i
and the values µ
i
where the function
is evaluated are
w
i
µ
i
0.555 555 55 −0.774 596 67
0.888 888 88 0.000 000 00
0.555 555 55 0.774 596 67
(Four digits behind the zero are enough for your work). To bring the integral
to the limits (−1, 1) use the transformation x = 1/2
[
(b − a)µ + b + a
]
.
4.6: Show that the δ-scaled phase function as defined by (4.35) is normalized.
4.7: Determine the system of algebraic equations that permit us to find Z (µ
i
),
see (4.54).
4.8: Check in detail the validity of (4.59) by substituting (4.56) into the relations
(4.57) and (4.58).
4.9: Consider the radiance expansion in the form (4.63). For the determination of
flux densities and for divergences only the term m = 0 is required. Substitute
(4.63) with m = 0 into (2.76) to obtain a set of differential equations in terms
of I
m=0
(τ ). Make full use of the orthogonality properties of the spherical
functions. You may wish to use the recursion formula
µP
l
(µ) =
l
2l + 1
P
l−1
(µ) +
l + 1
2l + 1
P
l+1
(µ)
132 Quasi-exact solution methods for the RTE
4.10: Well equipped with the experience of the previous problem, verify equation
(4.65). You may wish to employ the recursion formula
µP
m
l
(µ) =
l + m
2l + 1
P
m
l−1
(µ) +
l − m + 1
2l + 1
P
m
l+1
(µ)
4.11: Derive equation (4.92) from (4.89) by using (4.91) and other required defini-
tions. Only by carrying out the required derivation you will become familiar
with the notation of the finite difference method.
4.12: Analogously to the derivation of (4.97) find (4.98).
4.13: Use simple Cartesian vector operations to verify (4.130) and (4.131).
4.14: Verify (4.137) by carrying out the required comparison.
5
Radiative perturbation theory
5.1 Adjoint formulation of the radiative transfer equation
Problems in radiative transfer theory can be solved by various solution methods. In
the previous chapter we have discussed a number of these which we will classify as
the forward or the regular methods. In addition to applying the forward solutions,
it is also possible to use the so-called adjoint solution techniques which offer the
decisive advantage that for certain types of transfer problems the numerical effort
can be drastically reduced.
In this section we will formulate the adjoint technique. It will be necessary to
introduce a new terminology involving expressions such as the radiative effect and
the atmospheric response due to the presence of energy sources. By means of an
important but simple example we will demonstrate the numerical advantage that the
adjoint technique offers in comparison to the forward formulation. In Section 5.2 we
will introduce the perturbation technique and show how to apply it to the forward
as well as to the adjoint formulation.
The adjoint method originated as a purely mathematical tool for the solution of
linear operator equations. Discussions on this subject can be found in textbooks on
principles of applied mathematics such as Courant and Hilbert (1953), Friedman
(1956) and Keener (1988). Before presenting the basic radiative transfer theory in
the adjoint form, we would like to point out forcefully that this formulation is not
merely another solution method to solve the RTE. It is a method of reformulating
the transfer problem for maximum computational efficiency.
In the sequel, it will be of advantage to use the height z as the vertical coordinate
instead of the optical thickness. We start our analysis with the RTE for the total
radiation in the form, cf. (2.27)
1
d
ds
I (s, Ω) =−k
ext
(s)I (s, Ω) +
k
sca
(s)
4π
4π
P(s, Ω
→ Ω)I (s, Ω
)d
+ J
e
(s, Ω)
(5.1)
1
For brevity, the index ν will again be omitted.
133
134 Radiative perturbation theory
where we have used Ω ·∇I = dI/ds. Assuming horizontal homogeneity of all
variables the geometric increment ds may be replaced by dz/µ yielding
µ
∂
∂z
I (z, Ω) =−k
ext
(z)I (z, Ω) +
k
sca
(z)
4π
4π
P(z, Ω
→ Ω)I (z, Ω
)d
+J
e
(z, Ω)
(5.2)
We could have used the total derivative dI/dz since the variables µ and ϕ are
treated as constants. For simplicity we ignore thermal radiation and polarization
effects and consider the Sun as the only radiation source. If z = z
t
stands for the top
of the atmosphere then the source function for true emission J
e
(z, Ω), introduced
in (2.26), may be written as
J
e
(z, Ω) =|µ
0
|S
0
δ(z − z
t
)δ(µ − µ
0
)δ(ϕ − ϕ
0
) = Q(z, Ω) (5.3)
where δ denotes the Dirac δ-function. Equation (5.3) is by no means the only
possible formulation of the solar source. Had we used the RTE after the direct–
diffuse splitting of the radiation field a different form for the source function Q(z, Ω)
would have resulted.
We are now going to introduce the linear differential operator L
L = µ
∂
∂z
+ k
ext
(z) −
k
sca
(z)
4π
4π
P(z, Ω
→ Ω) ◦ d
(5.4)
which is a part of (5.2). The symbol ◦ occurring in the integral of this equation must
be replaced by the particular function on which the operator L operates. Utilizing
(5.3) and (5.4) the RTE can be written in the brief and elegant form
LI(z, Ω) = Q(z, Ω)
(5.5)
Before we discuss the so-called adjoint operator L
+
which is associated with L,
we need to introduce the definition of the inner product
f
1
, f
2
=
z
t
0
2π
0
1
−1
f
1
(z,µ,ϕ) f
2
(z,µ,ϕ)dµ dϕ dz
=
z
t
0
4π
f
1
(z, Ω) f
2
(z, Ω)d dz
(5.6)
where the bracket expression on the left-hand side is a shorthand notation for
the integral expressions of the right-hand side. Here f
1
and f
2
are two arbitrary
functions. We wish to emphasize that in the definition of the inner product the
integrals extend over the complete range of each variable, i.e. the integration extends
over the so-called phase space of a particular problem.