Zdunkowski W., Trautmann T., Bott A. Radiation in the atmosphere: A course in Theoretical Meteorology

Подождите немного. Документ загружается.

4.4 The spherical harmonics method 105

with m = 0, 1,..., and l = m, m + 1,...,M. In order to have as many

unknowns I

as equations we require

m−1

(τ ) = 0, I

M+1

(τ ) = 0 (4.66)

These two terms appear in (4.65) for l = m and l = M, respectively. Otherwise the

system would be under-determined.

We return to the discussion why there should be an even number (2 p) of equations

for the unknown functions I

(τ ). From (4.63) we see that the expansion coefﬁcients

(τ ) fully determine the m-th azimuthal Fourier expansion coefﬁcient I

(τ,µ)

of the radiance function. Considering an inhomogeneous atmosphere consisting of

individual homogeneous sublayers, for each m the radiance function I

(τ,µ)is

required to be a continuous function of µ in the range (−1, 1). Since at the top of a

particular layer we have the same number of determining equations for the incoming

downwelling radiance as for the incoming upwelling radiance at the bottom of that

layer, in total this leads to an even number of equations over the full range of µ.

For further details see also Dave (1975).

It is convenient to introduce the following vector and matrix symbols

(τ ) =







(τ )

m+1

(τ )







(4.67)

for m = 0, 1,...,. Note that this column vector has exactly 2p rows. For the

primary scattered sunlight and the thermal emission term, i.e. the terms on the

right-hand side of (4.65), we introduce the 2 p-dimensional vector f

(τ ). Further-

more, we need to deﬁne two 2p × 2 p coefﬁcient matrices A

, B

. The ﬁrst of

these contains the factors in (4.65) multiplying the derivatives with respect to

τ , while the second matrix incorporates the factors premultiplying I

. In this

manner we obtain from (4.65) for each Fourier mode m the matrix differential

equation

(τ )

dτ

− A

(τ ) = f

(τ ) (4.68)

The matrix elements of A

and B

can be determined as follows: let us consider a

ﬁxed but arbitrary value of m. The ﬁrst two and the last two equations of the system

106 Quasi-exact solution methods for the RTE

(4.65) are given by

(2m + 1)

m+1

(τ )

dτ

[

− (2m + 1)

]

(τ )

=−

2π

exp



−



(−µ

) − 2(1 − ω

)B(τ )δ

× (2m + 2)

m+2

(τ )

dτ

(τ )

dτ

[

m+1

− (2m + 3)

]

m+1

(τ )

=−

2π

exp



−



m+1

(−µ

) − 2(1 − ω

)B(τ )δ

0,m+1

(M + m)

(τ )

dτ

+ (M − 1 − m)

M−2

(τ )

dτ

[

M−1

− (2M − 1)

]

M−1

(τ )

=−

2π

exp



−



M−1

(−µ

) − 2(1 − ω

)B(τ )δ

0,M−1

× (M − m)

M−1

(τ )

dτ

[

− (2M + 1)

]

(τ )

=−

2π

exp



−



(−µ

) − 2(1 − ω

)B(τ )δ

0,M

(4.69)

It can be readily seen that the matrix B

is tridiagonal with zero entries on the

main diagonal. The matrix elements B

i,k

are given by

i,i−1

= i − 1, i = 2,...,2 p

i,i+1

= 2m + i, i = 1,...,2p − 1

i,k

= 0

i = 1,...,2 p

k = i − 1, i + 1

(4.70)

Similarly we ﬁnd for the matrix A

nonzero entries on the main diagonal

i,k

= 0, i = k,

i = 1,...,2 p

k = 1,...,2 p

i,i

= [2(m + i) −1] − ω

m+i−1

, i = 1,...,2 p

(4.71)

For the components of the vector f

(τ ) one obtains

(τ ) =−

2π

exp



−



i+m−1

(−µ

) − 2(1 − ω

)B(τ )δ

0,i+m−1

i = 1,...,2p (4.72)

4.4 The spherical harmonics method 107

Multiplication of (4.68) by (B

)

−1

and using the abbreviations

= (B

)

−1

, D

(τ ) = (B

)

−1

(τ ) (4.73)

yields

(τ )

dτ

− G

(τ ) = D

(τ ) (4.74)

The solution of this ordinary differential equation is

(τ ) = exp







exp[G

(τ − τ



)]D

(τ



)dτ



(4.75)

where C

is the vector of integration constants













(4.76)

It should be noted that the apparently complicated solution (4.75) contains parts

of the boundary conditions. For a homogeneous layer extending between τ = 0

and τ = τ

we directly obtain at the upper boundary

(τ = 0) = C

(4.77)

which determines the constants of integration after speciﬁcation of the radiation

incident at τ = 0. For more details the reader may consult the work of Flatau and

Stephens (1988).

There are many ways to evaluate the exponential matrix appearing in (4.75), see

Moler and van Loan (1978). One way is to determine the exponential matrix with

the help of the Jordan matrix J

in normal form which is deﬁned as







... 0

0 ... λ







(4.78)

The exponential matrix can be found from

exp(G

τ ) = P

exp(J

τ )(P

)

−1

(4.79)

where P

is the so-called modal matrix containing the eigenvectors of G

. The

individual eigenvectors result from the distinct eigenvalues λ

,...,λ

. A simple

discussion on the subject may be found in Bronson (1972).

108 Quasi-exact solution methods for the RTE

However, for systems not too large the Putzer method (Putzer, 1966) provides an

entirely analytical procedure which will be outlined next. The exponential matrix

occurring in (4.75) may be evaluated using the following algorithm

(τ ) = exp(G

τ ) =

M−m



j=0

j+1

(τ )X

(4.80)

In the previous formula we have used the deﬁnitions

= E, X

i=1



− λ



, j = 1, 2,...,M + 1 − m (4.81)

where E is the 2 p × 2 p unit matrix. The λ

are the eigenvalues of the matrix G

which can be determined using standard numerical routines. The scalar coefﬁcients

(τ ),...,η

M+1−m

(τ ) occurring in (4.80) can be recursively determined via the

solution of the following linear system of ordinary differential equations

dη

dτ

= λ

dη

j+1

dτ

= λ

j+1

+ η

, j = 1,...,M − m

(4.82)

The initial conditions for this system are given by

(0) = 1, η

j+1

(0) = 0 (4.83)

Although this method requires solving a system of linear differential equations, it

has a triangular structure. Thus the solutions can be determined in succession.

Using the deﬁnition (4.80) we may formally write for the integral term in (4.75)

(τ ) =



exp [G

(τ − τ



)]D

(τ



)dτ





(τ − τ



(τ



)dτ



(4.84)

Utilizing (4.80) and (4.84) the solution (4.75) of the m-th differential equation

system may be written as

(τ ) = N

(τ )C

+ H

(τ ) (4.85)

Alternatively one may write this equation in component form as







(τ )













1,1

(τ ) ··· N

1,2p

(τ )

2p,1

(τ ) ··· N

2p,2p

(τ )

























(τ )







(4.86)

From this equation we immediately see that for m = 0, 1,..., we have a total

of  + 1 different solutions. For each particular value of m there are in total 2 p

4.4 The spherical harmonics method 109

different rows in (4.86). It is noteworthy that for an increasing Fourier mode m the

system of differential equations decreases in size up to the point where for m = M

only a single scalar differential equation has to be solved.

We now return to (4.75) to involve the boundary conditions. The integration

constants C

are determined by ﬁrst combining the solutions I

(τ ) as contained in

(4.63) to obtain the m-th Fourier mode of the radiance. There is no prescribed way to

specify the boundary conditions. A method which has found wide acceptance is the

so-called Marshak boundary condition which will be used in the following. We will

assume that no downwelling diffuse radiation enters at the top of the atmosphere.

For the downwelling and the upwelling radiance ﬁelds I

(τ = 0, −µ) and I (τ

,µ)

the following relations must hold

(a)



−1

(τ = 0,µ)P

m+2 j −1

(µ) dµ = 0

(b)



(τ

,µ)P

m+2 j −1

(µ) dµ







−1

(τ

,µ



)µ



dµ





exp



−



+(1 − A





m+2 j −1

(µ)dµ (4.87)

where j = 1, 2,..., p and m = 0, 1,...,. Note that the scalar version of the

boundary condition (4.30) has already been employed in (4.87b). The integration

constants may be determined using standard numerical algorithms for systems of

linear equations. A vertically inhomogeneous atmosphere is handled as in DOM

by requiring continuity conditions for the radiation ﬁeld at the interior boundaries

for each layer.

The main advantages and disadvantages of the SHM are listed here.

(1) The solution of the radiance ﬁeld, to a large degree, can be derived in an analytic

manner.

(2) In contrast to DOM, a discretization of the µ-dependence is not required. Therefore,

the computation time does not increase when the radiance is needed for a large number

of directions µ.

(3) The SHM circumvents the problems involved in integrating the highly oscillatory P

functions for large l by carrying out these integrations analytically.

(4) The radiances at all depths inside as well as the reﬂected and transmitted radiation ﬁeld

may be obtained simultaneously.

(5) An increasing number of directions does not notably change the total computation time.

The same fact applies to the total optical depth of the medium.

(6) The SHM needs less computation time as MOM and SOS if the same amount of

information is required (internal radiation ﬁeld).

110 Quasi-exact solution methods for the RTE

(7) The SHM is closely related to DOM. Benchmark computations of both methods essen-

tially yield identical results.

(8) As is the case for any exact treatment of the RTE, the computation time drastically

increases with the number of expansion terms () in the phase function.

The SHM as formulated above follows the work of Deuze et al. (1973),

Zdunkowski and Korb (1974), and Zdunkowski and Korb (1985). It should also

be emphasized that for a total of four streams and for the azimuthally independent

radiation ﬁeld (m = 0) very efﬁcient and accurate solutions may be obtained for

ﬂuxes and heating rates, since both the eigenvalue problem as well as the vector

m=0

(τ ) can be solved analytically. For more details interested readers are directed

to the work of Li and Ramaswamy (1996) and Zdunkowski et al. (1998).

4.5 The ﬁnite difference method

The ﬁnite difference method (FDM) is based on the integro-differential form of the

RTE (2.76), using a phase function truncation after  terms, which for the m-th

Fourier mode is repeated below

dτ

(τ,µ) = I

(τ,µ) −



−1





l=m

(µ)P

(µ



(τ,µ



)dµ



−

4π

exp



−







l=m

(µ)P

(−µ

)−(1 − ω

)B(τ )δ

(4.88)

In the following it is more convenient to formulate the RTE in z-space rather than

the usual τ -space, i.e.

(z,µ) =−k

ext

(z)I

(z,µ)+

sca

(z)



−1





l=m

(µ)P

(µ



(z,µ



)dµ



sca

(z)

4π

exp



−

τ (z)







l=m

(µ)P

(−µ

) + k

abs

(z)B(z)δ

(4.89)

Let us discretize µ by introducing a discrete set of 2s Gaussian quadrature points

,µ

) with the usual properties

−i

= w

, µ

−i

=−µ

, i =−s,...,−1, 1,...,s (4.90)

4.5 The ﬁnite difference method 111

Furthermore, symmetric and antisymmetric sums of the radiation ﬁeld in direction

are introduced via

m,+

(z) =



(z) + I

−i

(z)



, I

m,−

(z) =



(z) − I

−i

(z)



(4.91)

where I

(z) = I

(z,µ

) and I

−i

(z) = I

(z, −µ

). Based on these deﬁnitions, and

after some tedious but simple steps, we obtain a set of coupled, ordinary ﬁrst-order

differential equations approximating the original transfer equation

(a) µ

m,+

+ k

ext

m,−

− k

sca



j=1

m,−

= k

sca

m,−

i = 1,...,s

(b) µ

m,−

+ k

ext

m,+

− k

sca



j=1

m,+

= k

sca

m,+

+ k

abs

Bδ

(4.92)

Here, S

m,+

, S

m,−

symbolize the primary scattered light. Note that for the sym-

metric sum in (4.92a) no Planckian emission term occurs. For the symmetric and

the antisymmetric sum of the phase function in directions ±µ

we have deﬁned

m,+

= p

(µ

,µ

) + p

(µ

, −µ

), P

m,−

= p

(µ

,µ

) − p

(µ

, −µ

)

(4.93)

where the coefﬁcients p

are given by

(µ

,µ

) =





l=m

(µ

) (4.94)

The deﬁnitions (4.93) make use of the symmetry properties of the phase function

P(−µ

, −µ

) = P(µ

,µ

), P(−µ

,µ

) = P(µ

, −µ

)

(4.95)

Finally, for the primary scattered sunlight the following expressions have been

introduced in (4.92)

m,+

4π

exp



−

τ (z)







l=m



(µ

) + P

(−µ

)



(−µ

)

m,−

4π

exp



−

τ (z)







l=m



(µ

) − P

(−µ

)



(−µ

)

(4.96)

The special case m = 0 is of particular interest if ﬂuxes and heating rates are

required. This is the situation treated by Barkstrom (1976) which will be discussed

next.

Since P

(−µ) = (−1)

(µ), see (2.59a), the Legendre expansion of the

The case m = 0 has also been investigated, see Rozanov et al. (1997).

112 Quasi-exact solution methods for the RTE

phase function allows us to write P

= P

m=0,+

and P

−

= P

m=0,−

as fully

symmetric expressions. We recognize this from the following expression.





l=0

(µ

) +





l=0

(µ

(−µ

)

/2



l=0

(µ

) +







l=0

2l+1

(µ

2l+1

(µ

)

/2



l=0

(µ

(−µ

) +







l=0

2l+1

(µ

2l+1

(−µ

)

/2



l=0

(µ

)

(4.97)

where 



is the next smallest integer value to ( − 1)/2 and p

= p

m=0

. Analo-

gously we obtain for P

−

the expression

−

/2



l=0

2l+1

(µ

2l+1

(µ

) (4.98)

In a similar manner we can deﬁne for the primary scattered sunlight the expressions

/2



l=0

(µ

), P

−

/2



l=0

2l+1

(µ

2l+1

(µ

)

(4.99)

so that the terms involving the primary scattered direct sunlight are given by

4π

exp



−

τ (z)



, S

−

4π

exp



−

τ (z)



−

(4.100)

4.5.1 Vertical discretization in the ﬁnite difference method

In order to solve the coupled system of differential equations (4.92) for m = 0, we

introduce an odd number of discrete values z

, k = 1, 2,...,2K + 1. The deriva-

tives occurring in (4.92) will be approximated by means of centered differences



≈

f (z

k+1

) − f (z

k−1

)

z

with z

= z

k+1

− z

k−1

(4.101)

For numerical convenience Barkstrom (1976) recommends to evaluate I

−

m=0,−

at all even-numbered grid points and at the boundaries z

, z

2K +1

, and

4.5 The ﬁnite difference method 113

the function I

= I

m=0,+

at all odd-numbered grid points. Experience shows that

these approximations produce rather accurate results for the case of diffuse incident

radiation.

In case of an incident direct solar beam the primary scattered light varies very

rapidly with optical depth so that the above procedure is not sufﬁciently accurate.

However, this deﬁciency can be eliminated by analytically integrating the primary

scattered light as will be shown below.

Integrating (4.92) between the interior points z

k−1

and z

k+1

and approximating

all integrals (excepting the primary scattering term) by means of



k+1

k−1

f (z)dz ≈ z

f (z

) (4.102)

gives

k+1

) − I

k−1

) +

z



ext

−

) − k

sca

)



j=1

−

)





k+1

k−1

sca

(z)

−

(z)dz

−

k+1

) − I

−

k−1

) +

z



ext

) − k

sca

)



j=1

)





k+1

k−1

sca

(z)

(z)dz +

z

abs

)B(z

) (4.103)

Due to the identity

exp



−

τ (z)



ext

(z)

exp



−

τ (z)



(4.104)

we obtain



k+1

k−1

sca

(z)

(z)dz ≈

sca

)

ext

)

4π

)



exp



−

τ (z

k+1

)



− exp



−

τ (z

k−1

)



(4.105)

As already mentioned above, (4.105) provides a very accurate approximation for

the integration of the solar term.

114 Quasi-exact solution methods for the RTE

Substituting (4.105) into (4.103) yields the discretized form of the RTE

k+1

) − I

k−1

)

z



ext

−

) − k

sca

)



j=1

−

)



sca

)

ext

)

4π

−

)



exp



−

τ (z

k+1

)



− exp



−

τ (z

k−1

)



−

k+1

) − I

−

k−1

)

z



ext

) − k

sca

)



j=1

)



sca

)

ext

)

4π

)



exp



−

τ (z

k+1

)



− exp



−

τ (z

k−1

)



z

abs

)B(z

)

(4.106)

4.5.2 Treatment of the boundary conditions

Next we will consider the speciﬁcation of the boundary conditions. We will assume

that no diffuse radiation is incident at the top of the atmosphere, i.e. at z = z

2K +1

I (z

,µ < 0) = 0 (4.107)

At the lower boundary z = z

= 0 we assume an isotropically emitting ground

with temperature T

. However, in accordance with observations the ground emits

only a fraction ε

of the black body radiation B(T

). The term ε

is called the

emissivity of the ground. This concept will be discussed in more detail in a later

chapter. Consequently, we must allow for isotropic reﬂection of diffuse radiation

with albedo A

. The boundary condition then is given by

I (0,µ > 0) = 2A



I (0, −µ



)µ



dµ



exp



−

τ (0)



+ ε

B(T

)

(4.108)

According to (4.91) the upwelling and downwelling diffuse radiation can be

recovered from the I

- and I

−

-functions via

(z) = I

(z) + I

−

(z), I

−i

(z) = I

(z) − I

−

(z), i = 1,...,s (4.109)