Zdunkowski W., Trautmann T., Bott A. Radiation in the atmosphere: A course in Theoretical Meteorology
Подождите немного. Документ загружается.
5.5 Appendix 155
which includes the total downward flux density E
v,−
(τ
0
) resulting from the vacuum
problem (i), that is
E
v,−
(τ
0
) =
2π
0
0
−1
|µ
|I
v
(τ
0
,µ
,ϕ
)dµ
dϕ
+|µ
0
|S
0
exp
−
τ
0
|µ
0
|
(5.94)
Note that in (5.92) the factor 1/π converts the flux density to the isotropic radiance.
Let us now consider another special radiative transfer problem in terms of a
dimensionless radiance I
g
assuming an isotropic radiation source of unit strength at
the lower boundary and vacuum boundary condition at the top of the atmosphere.
This transfer problem is defined by
µ
d
dτ
I
g
(τ,µ) = I
g
(τ,µ) −
ω
0
2
1
−1
P(τ, µ
,µ)I
g
(τ,µ
)dµ
(5.95)
with the boundary conditions
top of the atmosphere: I
g
(0,µ) = 0, −1 ≤ µ<0
Earth’s surface: I
g
(τ
0
,µ) = 1, 0 <µ≤ 1
(5.96)
The solution to this problem will be employed to find the solution to the RTE (5.91).
This is possible since I
g
(τ,µ) represents the response of the atmosphere due to an
isotropic diffuse illumination at τ = τ
0
with flux density
E
g,+
(τ
0
) = 2π
1
0
µI
g
(τ
0
,µ)dµ = π (5.97)
The total flux density in upward direction of the combined problem (i) and (ii)
is given by
E
+
(τ
0
) = A
g
E
−
(τ
0
) (5.98)
so that
E
+
(τ
0
)
E
g,+
(τ
0
)
=
A
g
π
E
−
(τ
0
) (5.99)
Analogously to this expression for I
r
d
(τ,µ) we may write
I
r
d
(τ,µ)
I
g
(τ,µ)
=
A
g
π
E
−
(τ
0
)orI
r
d
(τ,µ) =
A
g
π
E
−
(τ
0
)I
g
(τ,µ) (5.100)
In order to comprehend this conclusion, first of all, we recognize the similarity
of the transfer equations (5.91) and (5.95) and the formal agreement between the
upper boundary conditions in (5.92) and (5.96). The lower boundary conditions in
(5.92) and (5.96) vary in form, but both refer to isotropic upward radiation. The
step leading to (5.100) is motivated by the fact that the transfer problem for I
g
is
linear with respect to the isotropic boundary conditions as applied to I
g
at τ = τ
0
.
156 Radiative perturbation theory
This implies that a multiple of the upwelling radiation illuminating the atmosphere
from below results in a corresponding multiple of I
g
(τ,µ) which represents the
response of the atmosphere to the lower boundary conditions in (5.92).
To complete the problem, we need to express the factor multiplying I
g
(τ,µ)in
(5.100) in terms of known quantities. With the help of (5.94) and (5.100) we obtain
A
g
π
E
−
(τ
0
) =
A
g
π
E
v,−
(τ
0
) + 2π
0
−1
|µ
|I
r
d
(τ
0
,µ
)dµ
=
A
g
π
E
v,−
(τ
0
) + 2A
g
E
−
(τ
0
)
0
−1
|µ
|I
g
(τ
0
,µ
)dµ
(5.101)
By introducing the abbreviation
E
g,−
(τ
0
) = 2π
0
−1
|µ
|I
g
(τ
0
,µ
)dµ
(5.102)
into (5.101) we find the required expression
E
−
(τ
0
) =
1
1 − (A
g
/π)E
g,−
(τ
0
)
E
v,−
(τ
0
)
= E
v,−
(τ
0
)
1 +
A
g
π
E
g,−
(τ
0
) +
A
g
π
E
g,−
(τ
0
)
2
+···
(5.103)
The physical interpretation of this expression is as follows: the first term represents
the downward flux density of problem (i), while terms two, three, etc., describe the
additional contributions due to one, two, etc., reflection interactions between the
Lambertian surface and the plane–parallel atmosphere.
Substituting (5.103) into (5.100) gives
I
r
d
(τ,µ) = E
v,−
(τ
0
)
(A
g
/π)
1 − (A
g
/π)E
g,−
(τ
0
)
I
g
(τ,µ) (5.104)
The radiance I (τ,µ,ϕ) is the combination of the radiances I
v
(τ,µ,ϕ) and I
r
d
(τ,µ),
i.e.
I (τ,µ,ϕ) = I
v
(τ,µ,ϕ) + I
r
d
(τ,µ)
= I
v
(τ,µ,ϕ) + E
v,−
(τ
0
)
(A
g
/π)
1 − (A
g
/π)E
g,−
(τ
0
)
I
g
(τ,µ)
(5.105)
While in (5.32) we have written the general form for I
g
(z,µ,ϕ), here we have
omitted the azimuthal dependency of this quantity which is the actual form used
by Box et al. (1989b) in their computations.
Finally, it is noteworthy that in Chapter 4 we have included an equation similar to
(5.32), see (3.48), which is due to Chandrasekhar (1960). We wish to make reference
5.6 Problems 157
to Muldashev et al. (1999) who have generalized Chandrasekhar’s approach by
deriving the internal radiance field at an arbitrary level τ due to the presence of a
Lambertian surface as described in detail above.
5.6 Problems
5.1: Show that the third terms in (5.16) cancel out.
5.2: Write down in detail all steps involved in (5.64).
5.3: Verify both parts in (5.72).
5.4: Show that (5.73) is correct.
6
Two-stream methods for the solution
of the radiative transfer equation
6.1 δ-scaling of the phase function
A particular difficulty when solving the RTE for solar radiation in an atmosphere
containing aerosol particles and cloud droplets is associated with the fact that the
scattering phase function for such particle populations is highly peaked in the
forward direction. For cloud particles the energy scattered within an angle interval
of about 5
◦
around the forward direction is four to five magnitudes larger than that
part of the energy related to sideward or backward scattering. Therefore, it can
be concluded that an accurate treatment of the radiation field for highly peaked
phase functions requires a very large number of expansion terms when, according
to (2.55), the phase function is written as a series of Legendre polynomials.
In order to incorporate the forward diffraction contribution in multiple scattering
computations we may proceed as follows. Let f represent that part of the radiation
which is scattered in the forward direction. In specifying this forward part we may
define the so-called δ-scaled phase function by
P
∗
(cos ) = 2 f δ(1 − cos ) + (1 − f )
n−1
k=0
p
∗
k
P
k
(cos ) (6.1)
where the first term on the right-hand side represents the forward scattering
contribution. The modified expansion coefficients p
∗
k
will now be determined in
such a way that the first n moments p
∗
0
, p
∗
1
,..., p
∗
n−1
of the modified phase func-
tion P
∗
are equal to the corresponding moments of the original phase function P.
With the help of (2.61) we may, therefore, write
1
−1
P
l
(cos )P
∗
(cos )d cos =
1
−1
P
l
(cos )P(cos )d cos , l =0,...,n − 1
(6.2)
158
6.1 δ-scaling of the phase function 159
According to Morse and Feshbach (1953) the δ-function can be expressed as an
infinite series of Legendre polynomials, that is
δ(x − x
) =
∞
k=0
2k + 1
2
P
k
(x)P
k
(x
) (6.3)
Utilizing this expression together with the orthogonality relations for the Legendre
polynomials (2.59b), we obtain
1
−1
P
l
(cos )2 f δ(1 − cos )d cos
=
1
−1
P
l
(cos )2 f
∞
k=0
2k + 1
2
P
k
(1)P
k
(cos )d cos = 2 f (6.4)
since P
k
(1) = 1, see Abramowitz and Stegun (1972). Hence, (6.2) may be rewritten
as
2 f + (1 − f )
1
−1
P
l
(cos )
n−1
k=0
p
∗
k
P
k
(cos )d cos
=
1
−1
P
l
(cos )P(cos )d cos l = 0, 1,...,n − 1 (6.5)
By applying once more the orthogonality relations (2.59b) we finally obtain
p
∗
l
=
p
l
− (2l + 1) f
1 − f
, l = 0, 1,...,n − 1
(6.6)
So far we have given no rule for determining f . Since the modified phase function
is truncated after the (n − 1)th term, that is p
∗
l
= 0 for l ≥ n, we obtain from (6.6)
a relation between f and the p
l
f =
p
l
2l + 1
, l ≥ n (6.7)
However, this expression does not determine f uniquely. In order to come up with
a suitable choice we consider the difference between the original phase function P
and the scaled phase function P
∗
. In noting that the δ-function can be expanded in
an infinite series of Legendre polynomials, see (6.3), the difference between P and
P
∗
can be determined from
P(cos ) − P
∗
(cos ) =
∞
l=n
[
p
l
− (2l + 1) f
]
P
l
(cos ) (6.8)
By choosing
f =
p
n
2n + 1
(6.9)
160 Two-stream methods for the solution of the RTE
we observe that the leading term in the series (6.8) vanishes, while the remaining
terms of the series do not vanish. On the basis of numerical experience, Wiscombe
(1977) argues that it is more important to force agreement in the phase function for
the lowest order term l = n than for higher order terms l > n. Therefore, we will
use (6.9) to determine the fraction of the radiation which is scattered in the forward
direction.
In the following we will demonstrate how the δ-scaling of the phase function
affects the RTE in the absence of thermal emission. With reference to (2.29) and
(2.47) we write the RTE for the total radiance I
tot
including the parallel solar
radiation as
dI
tot
ds
=−k
ext
I
tot
+ k
ext
J (6.10)
where the multiple scattering term is given by
J =
ω
0
4π
4π
P(cos )I
tot
(s, Ω
)d
(6.11)
Introducing the δ-approximation of the phase function (6.1) into the multiple scat-
tering term and noting that
δ(1 − cos ) = 2πδ(µ − µ
)δ(ϕ − ϕ
) (6.12)
(see, e.g. Jackson, 1975), we obtain
J =
ω
0
4π
2π
0
1
−1
4π f δ(µ − µ
)δ(ϕ − ϕ
) + (1 − f )
n−1
l=0
p
∗
l
P
l
(cos )
× I
tot
(s,µ
,ϕ
)dµ
dϕ
(6.13)
The term containing the two δ-functions can be evaluated immediately so that
J = ω
0
fI
tot
(s,µ,ϕ) +
ω
0
(1 − f )
4π
2π
0
1
−1
n−1
l=0
p
∗
l
P
l
(cos )I
tot
(s,µ
,ϕ
)dµ
dϕ
(6.14)
For the RTE we then obtain
dI
tot
ds
=−k
∗
ext
I
tot
+ k
∗
ext
J
∗
(6.15)
where the modified source term J
∗
is given by
J
∗
=
ω
∗
0
4π
2π
0
1
−1
n−1
l=0
p
∗
l
P
l
(cos )I
tot
(s,µ
,ϕ
)dµ
dϕ
(6.16)
The modified extinction coefficient k
∗
ext
and the modified single scattering albedo
ω
∗
0
are
k
∗
ext
= (1 − ω
0
f )k
ext
, ω
∗
0
=
(1 − f )ω
0
1 − ω
0
f
(6.17)
6.2 The two-stream radiative transfer equation 161
By comparing (6.15) with (6.10) we observe that both types of the radiative
transfer equation have the same form. This means that the δ-adjustment of the
phase function leaves the RTE formally invariant. In summary, in going from the
original unscaled to the scaled problem we have to make the following replacements
p
l
→ p
∗
l
=
p
l
− (2l + 1) f
1 − f
k
ext
→ k
∗
ext
= (1 − ω
0
f )k
ext
ω
0
→ ω
∗
0
=
(1 − f )ω
0
1 − ω
0
f
(6.18)
where the forward diffraction part f of the phase function is given by (6.9).
For a homogeneous layer with thickness z the optical depth τ = k
ext
z has
to be replaced by the scaled value
τ
∗
= k
∗
ext
z = (1 − ω
0
f )τ
(6.19)
It is important to note that in the δ-approximation for consistency the extinction of
the direct solar beam is given by
S(τ
∗
) = S
0
exp
−
τ
∗
µ
0
(6.20)
Since τ
∗
≤ τ , the δ-scaled direct flux density is always larger than its unscaled
counterpart. This means that the sum of the scaled direct beam plus the scaled diffuse
radiation should be compared with the measured values of the global radiation, that
is the total downward radiation.
The δ-scaling approximation for solar radiative transfer has been widely
employed. In particular, this method has been applied to the so-called two-stream
methods, see Joseph et al. (1976) and Wiscombe (1977). These authors have shown
that the δ-scaling resulted in improved flux density and heating rate calculations.
6.2 The two-stream radiative transfer equation
It is well known that in a plane–parallel atmosphere the evaluation of the thermo-
dynamic heat equation requires the knowledge of the net radiative flux densities
E
net,z
= E
+,z
− E
−,z
but does not need the complete directional dependence of the
radiation field.
1
An accurate determination of these flux densities requires an inte-
gration of the azimuthally independent radiation field I
m=0
(τ,µ) over all directions
as shown in (2.126).
In many circumstances the computational effort to determine the up- and
downward flux densities in this manner is far too high so that approximate methods
1
For simplicity, the subscript z occurring at the radiative flux densities will henceforth be omitted.
162 Two-stream methods for the solution of the RTE
must be used. These are the so-called two-stream methods (TSM) which originated
with Schuster (1905), Schwarzschild (1906), Emden (1913), and Eddington (1916).
Owing to their high computational efficiency, two-stream methods are now widely
used in conjunction with global or mesoscale climate models for determining the
heating and cooling rates in the solar and long-wave radiative spectrum. In the fol-
lowing we will show how various versions of the TSM can be derived from the RTE.
In principle we could start from DOM by considering two radiation streams only.
However, it is more expedient to begin with the azimuthally averaged radiation field
in the absence of thermal radiation since the solar and infrared spectra may easily be
separated. For convenience we repeat (2.81), omitting there the infrared radiation
term
µ
d
dτ
I (τ,µ) = I (τ, µ) −
ω
0
2
1
−1
P(µ, µ
)I (τ,µ
)dµ
−
ω
0
4π
S
0
exp
−
τ
µ
0
P(µ, −µ
0
) (6.21)
The up- and downward radiative flux densities follow from
E
+
= 2π
1
0
µI (τ,µ)dµ, E
−
= 2π
1
0
µI (τ,−µ)dµ (6.22)
so that the total net flux density is obtained from
E
net
= E
+
− E
−
− µ
0
S
0
exp
−
τ
µ
0
(6.23)
see also Section 2.7. The last term on the right-hand side of (6.23) represents the
contribution of the unscattered solar beam. Integration of (6.21) over the 2π upper
hemisphere leads to
dE
+
dτ
= 2π
1
0
I (τ,µ)dµ − ω
0
π
1
0
1
−1
P(µ, µ
)I (τ,µ
)dµ
dµ
− S
0
exp
−
τ
µ
0
ω
0
2
1
0
P(µ, −µ
0
)dµ (6.24)
Let us now discuss that part of the diffuse downward light which is scattered
into the backward hemisphere. From the phase function we define the so-called
backscattered fraction or backscattering coefficient b(−µ
) (Wiscombe and Grams,
1976) as
b(−µ
) =
1
2
1
0
P(µ, −µ
)dµ, µ
> 0 (6.25)
where, as usual, −µ
refers to the direction of the downward diffuse light as
illustrated in Figure 6.1. Due to the normalization condition (4.34) and the symmetry
6.2 The two-stream radiative transfer equation 163
Upper hemisphere
Lower hemisphere
I(
τ, −µ
′
)
b(−µ
′
)
1 − b(−µ
′
)
Fig. 6.1 Illustration of the backscattered fraction b(−µ
) in the phase function.
properties (4.95) of P(µ, µ
) it follows that
1
2
1
0
P(−µ, −µ
)dµ =
1
2
1
0
P(µ, µ
)dµ = 1 − b(−µ
), µ
> 0 (6.26)
With this definition we can split the multiple scattering term in (6.24) into two parts
1
2
1
0
1
−1
P(µ, µ
)I (τ,µ
)dµ
dµ
=
1
2
1
0
1
0
P(µ, µ
)I (τ,µ
)dµ
dµ +
1
2
1
0
0
−1
P(µ, µ
)I (τ,µ
)dµ
dµ
=
1
2
1
0
1
0
P(µ, µ
)I (τ,µ
)dµ
dµ +
1
2
1
0
1
0
P(µ, −µ
)I (τ,−µ
)dµ
dµ
=
1
0
[1 − b(−µ
)]I (τ,µ
)dµ
+
1
0
b(−µ
)I (τ,−µ
)dµ
(6.27)
which motivates the definition (6.25). For the primary scattered sunlight in the
backward direction we obtain analogously to (6.25)
b(−µ
0
) =
1
2
1
0
P(µ, −µ
0
)dµ, µ
0
> 0 (6.28)
Substituting (6.25)–(6.28) into (6.24) yields
dE
+
dτ
= 2π
1
0
I (τ,µ)dµ − 2πω
0
1
0
[1 − b(−µ
)]I (τ,µ
)dµ
−2πω
0
1
0
b(−µ
)I (τ,−µ
)dµ
− ω
0
S
0
exp
−
τ
µ
0
b(−µ
0
) (6.29)
164 Two-stream methods for the solution of the RTE
In an analogous manner we can derive the corresponding equation for the downward
flux density E
−
dE
−
dτ
=−2π
1
0
I (τ,−µ)dµ + 2πω
0
1
0
b(−µ
)I (τ,µ
)dµ
+ 2πω
0
1
0
[1 − b(−µ
)]I (τ,−µ
)dµ
+ ω
0
S
0
exp
−
τ
µ
0
[
1 − b(−µ
0
)
]
(6.30)
Next we define the average of a quantity ψ(τ, µ) weighted with respect to the
up- and downwelling intensities I (τ,±µ)
ψ
+
(τ ) =
1
0
ψ(τ, µ)I (τ,µ)dµ
1
0
I (τ,µ)dµ
, ψ
−
(τ )=
1
0
ψ(τ, µ)I (τ,−µ)dµ
1
0
I (τ,−µ)dµ
(6.31)
For the particular choice ψ = µ we obtain from (6.31)
µ
+
(τ ) =
2π
1
0
µI (τ,µ)dµ
2π
1
0
I (τ,µ)dµ
=
E
+
2π
1
0
I (τ,µ)dµ
µ
−
(τ ) =
2π
1
0
µI (τ,−µ)dµ
2π
1
0
I (τ,−µ)dµ
=
E
−
2π
1
0
I (τ,−µ)dµ
(6.32)
The parameters µ
+
and µ
−
can be interpreted as the mean directional cosines of
the up- and downward diffuse radiation. Utilizing in (6.31) ψ = b(−µ) we obtain
in the same way
b
+
(τ ) =
2π
1
0
b(−µ)I (τ, µ)dµ
2π
1
0
I (τ,µ)dµ
, b
−
(τ ) =
2π
1
0
b(−µ)I (τ, −µ)dµ
2π
1
0
I (τ,−µ)dµ
(6.33)
It should be emphasized that, in general, the parameters µ
±
, b
±
are functions of
the optical depth τ since the radiances depend on τ . Combining (6.32) and (6.33)