A Monte Carlo method for 3D thermal infrared radiative transfer

ARTICLE IN PRESS

Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178 www.elsevier.com/locate/jqsrt

A Monte Carlo method for 3D thermal infrared radiative transfer Y. Chen, K.N. Liou Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles Los Angeles, CA 90095, USA Received 11 April 2005; received in revised form 27 July 2005; accepted 2 October 2005

Abstract A 3D Monte Carlo model for specific application to the broadband thermal radiative transfer has been developed in which the emissivities for gases and cloud particles are parameterized by using a single cubic element as the building block in 3D space. For spectral integration in the thermal infrared, the correlated k-distribution method has been used for the sorting of gaseous absorption lines in multiple-scattering atmospheres involving 3D clouds. To check the Monte-Carlo simulation, we compare a variety of 1D broadband atmospheric fluxes and heating rates to those computed from the conventional plane-parallel (PP) model and demonstrate excellent agreement between the two. Comparisons of the Monte Carlo results for broadband thermal cooling rates in 3D clouds to those computed from the delta-diffusion approximation for 3D radiative transfer and the independent pixel-by-pixel approximation are subsequently carried out to understand the relative merits of these approaches. r 2006 Elsevier Ltd. All rights reserved. Keywords: 3D thermal radiative transfer; Monte Carlo method; Cubic element

1. Introduction The Monte Carlo method has long been employed to simulate the transfer of radiation from the sun in planetary atmospheres [1]. Application to the problem involving the interaction of solar radiation with 3D cloud fields have also been made by many researchers [2–5]. The Monte Carlo method by virtue of its flexibility in tracing photons has been thought of as the exact solution for radiative transfer in multidimensional space. The Monte Carlo results have been frequently used as bench marks for other analytical 3D radiative transfer approaches [6–9]. Despite significant advances in the computer technology in terms of speed and storage, the ‘‘exact’’ Monte Carlo calculations are extremely expensive because a large number of photons must be employed to reduce the noise inherent in the use of randomly selected variables. A review of the literature reveals that the Monte Carlo approach has not been rigorously developed for broadband thermal infrared radiative transfer in atmospheres containing 3D clouds. In the field of neutron transport, emission from the hot interior is an important theoretical subject [10]. However, the photon Corresponding author. Tel.: +1 310 825 5039; fax: +1 310 267 2252.

E-mail address: [email protected] (K.N. Liou). 0022-4073/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2005.10.002

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

167

emission has been treated only as a function of temperature without accounting for the absorption and emission properties of pertinent gases. A Monte Carlo simulation of thermal infrared radiative transfer in 3D clouds has also been presented [11], but the emitted photons were assigned rather than physically formulated. Several papers have also reported the Monte Carlo simulations of microwave radiative transfer by means of a backward method combined with a forward method for radiance calculations at the satellite position [12–14]. We have developed a 3D broadband thermal infrared Monte Carlo radiative transfer model by combining a 3D Monte Carlo photon transport scheme with emission initiated from the medium, and the spectral radiation model developed by Fu and Liou [15,16], as modified in Chen et al. [9] for 3D radiative transfer. The specific features of this Monte Carlo model include the following: parameterization of the gaseous absorption and cloud emissivity for 3D radiative transfer; construction of the complex cloud geometry using the basic cubic cell approach; computation of the broadband atmospheric flux at each cell surface; and spectral integration using the correlated k-distribution method for the gaseous absorption in multiple-scattering atmospheres involving 3D clouds. This paper is organized as follows. In Section 2 we present the basic Monte Carlo method for thermal emission and the manner in which the gaseous absorption and cloud particles’ emissivity are parameterized. In Section 3 we compare the thermal infrared fluxes and cooling rates produced from the present model to those computed from the conventional plane-parallel approach. Additional comparisons have also been conducted with reference to the 3D delta-diffusion and pixel-by-pixel approaches. Summary is given in Section 4. 2. A Monte Carlo model for 3D thermal radiative transfer 2.1. Theory and model description The domain for Monte Carlo simulations in this model consists of a variable 3D structure composed of cubic grid cells, each of which is considered to represent a homogeneous and isothermal volume. All matterbased grid quantities such as the extinction coefficient bext, the single scattering albedo $, the phase function PðYÞ, the temperature T, and the emissivity e are stored locally in the cell center. A photon life in the thermal emission Monte Carlo program is essentially the same as in the solar counter part, except that the photons are emitted from grid cells and are not from the top of atmosphere (TOA). Thus, instead of looping over the emitting photons from the sun, looping over the grid cells has been performed in which the emitting photons are a function of emissivity and temperature. For a given wavelength l, the cubic grid cell (i,j,k) is assumed to have an emissivity el(i,j,k) and a temperature T(i,j,k), while the surface boundary cell (i,j) (two-dimension only) is defined by an emissivity esl(i,j) and a temperature Ts(i,j). The Planck functions for the cell (i,j,k) and the surface cell (i,j) are denoted by Bl(T(i,j,k)) and Bl(Ts(i,j)), respectively. The emitted radiance from each of the six cubic sides consisting of homogeneous and isothermal material is isotropic and can be written in the form I l ði; j; kÞ ¼ el ði; j; kÞBl ðTði; j; kÞÞ.

(1a)

The emitted radiance from the surface assuming Lambertian can also be written as I sl ði; jÞ ¼ esl ði; jÞBl ðT s ði; jÞÞ.

(1b)

The emitted fluxes from six cubic sides are defined by [17], see also the fluxes defined in Eqs. (16)–(21) Z Iði; j; k; OÞOxi dO F
xi ði; j; kÞ ¼ 2p

¼ pIði; j; kÞ,

ð2aÞ

where x1 ¼ x, x2 ¼ y, x3 ¼ z, Ox ¼ ð1  cos2 yÞ1=2 cos f, Oy ¼ ð1  cos2 yÞ1=2 sin f, and Oz ¼ cos y. Angular integrations over the upward and downward hemispheres are the same as in the PP case, i.e., (0, 2p) for f and (0,71) for cos y. Angular integrations in the x-direction are (p/2, p/2) and (1, 1), and (p/2, 3p/2) and (1, 1). In the y-direction, angular integrations are (0, p) and (1, 1), and (p, 2p) and (1, 1). The resulting integration is given in the second line of Eq. (2a).

ARTICLE IN PRESS 168

Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

The total flux for a given wavelength covering the model domain can then be expressed by XXX XX pI l ði; j; kÞ þ pI sl ði; jÞ, F total;l ¼ 6 i

j

k

i

(2b)

j

where 6 accounts for emission from the six cubic surfaces. If we assign a total photon number Ntotal to the model domain, we can relate it to the total flux integrated over the entire infrared spectrum, Ftotal. The number of photons for a cubic grid cell (i,j,k) and a surface cell (i,j) for a given wavelength l can be expressed by N l ði; j; kÞ ¼ 6 N sl ði; jÞ ¼

N total pel ði; j; kÞBl ðTði; j; kÞÞ, F total

N total pesl ði; jÞBl ðT s ði; jÞÞ. F total

(3)

(4)

The thermal emission from the surface boundary cell (i,j) is in the upward direction so that the photon’s emission direction is given by pffiffiffiffiffiffi (5) m ¼ cos y ¼  x1 ; f ¼ 2px2 , where x1 and x2 are random numbers between 0 and 1, and the negative value for m refers to the upward direction. The thermal emission from the cubic grid cell (i,j,k) is emitted from the six cell surfaces such that each emits Nl(i,j,k)/6 photons by first selecting a random location on the cell surface and then an isotropic emission direction. Each photon is assigned an initial weight of W m 0 ¼ 1, where m is the mth photon. Basically, the present Monte Carlo program includes three associative parts: the input atmospheric profile and cloud information that are used in calculation of the extinction coefficient bext, the single scattering albedo $, and the phase function PðYÞ for each cell (i,j,k); parameterization of the gaseous and cloud emissivities for individual cubic cells; and tracing each photon emitted from the cell surface until it is totally absorbed or leaves the model domain. The flux variable, F m n , which keeps track of the weight of the photon emitted from the cell surface, is defined by m Fm n¼0 ¼ W 0 ,

(6)

m where W m 0 defined previously is equal to 1 initially. After scattering and absorption events, W n (n ¼ 1; 2; . . .) decreases accordingly. The optical depth for the first scattering event is given by

tm 1 ¼  lnð1  xÞ,

(7)

where x is a random number sampled uniformly from 0 to 1. Note that ln (1x) has the same probability distribution as ln x. If tm 1 is greater than the cumulated optical depth the photon that travels from the starting point to the boundary surface, then whether the photon escapes the model domain or not depends on the boundary condition. If the photon remains in the domain, it is scattered in a direction determined by the cumulative probability function of the phase function at the scattering site. The weight of the photon after the first scattering is m Wm 1 ¼ $W 0 .

(8)

The fraction of the photon that is absorbed at the scattering site is m Am 1 ¼ ð1  $ÞW 0 .

(9)

The flux variable corresponding to the fraction of the photon that is scattered at first scattering site and crosses the site surface is  m W 1 if photon crosses the surface; Fm (10) ¼ 1 0 otherwise: The next scattering direction ðy; fÞ is determined by the phase function of the scattering site. The nth scattering quantities can also be obtained similar to Eqs. (7)–(10).

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

169

The optical depth to the nth scattering, tm n , is determined by tm n ¼  lnð1  xÞ,

(11)

which holds for nX1. The weight of the photon after the nth scattering is m Wm n ¼ $W n1 ,

(12)

which holds for nX1. The fraction of the photon that is absorbed at the scattering site is m Am n ¼ ð1  $ÞW n1 ,

(13)

which holds for nX1. The flux variable corresponding to the fraction of the photon that is scattered at the nth scattering site and crosses the site surface is given by  m W n if photon crosses the surface; Fm ¼ (14) n 0 otherwise: where we note that F m n is one of the six grid-cell surface flux variables. The total number of scatterings a photon undergoes is limited and is less than the prescribed nmax . The photon is terminated after the weight is less than the prescribed threshold value. We can convert the total weight of photons that are absorbed in each cell and across each cell surface into flux according to the relationship denoted in Eqs. (3) and (4). The 3D absorbed flux matrix is calculated from the following equation: Z F total l2 X X m Aði; j; kÞ ¼ A ði; j; kÞ dl, (15) N total l1 m n n;l where the sums over m and n are only done over the scatterings that occur in cell ði; j; kÞ and Am n is the energy that is absorbed at the nth scattering site of the mth photon at the wavelength l. The flux for each surface of cell (i,j,k) may be expressed by Z F total l2 X X "m " F z ði; j; kÞ ¼ F ði; j; kÞ dl, (16) N total l1 m n z;n;l F #z ði; j; kÞ ¼

F total N total

Z

F total F x ði; j; kÞ ¼ N total F! x ði; j; kÞ ¼

F y ði; j; kÞ ¼

F! y ði; j; kÞ

F total N total F total N total

F total ¼ N total

l2

XX

l1

Z

m l2

l1

Z

l2 l1

Z

l2 l1

Z

l2 l1

F #m z;n;l ði; j; kÞ dl,

XX m

F !m x;n;l ði; j; kÞ dl,

(19)

m F y;n;l ði; j; kÞ dl,

(20)

F !m y;n;l ði; j; kÞ dl,

(21)

n

XX m

(18)

n

XX m

m F x;n;l ði; j; kÞ dl,

n

XX m

(17)

n

n

where the sums over m and n are performed over the emission and scattering that occur in cell (i,j,k), and F "m z;n;l , m !m m !m " # ! ! F #m surfaces, z;n;l , F x;n;l , F x;n;l , F y;n;l , andF y;n;l are the flux variables that cross the z , z , x ,x , y , and y respectively. Note that n ¼ 0 for emission and nX1 for the scattering event of the mth photon for a given wavelength l.

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

170

2.2. Parameterization of gaseous and ice cloud emissivities An essential step in the development of the present 3D thermal radiative transfer model is parameterization of the grid cell emissivity. We have developed accurate parameterizations for the gaseous emissivity in terms of optical depth and the ice cloud emissivity in terms of absorption optical depth in PP condition and in a single 3D cubic cell, respectively. In the thermal infrared spectrum (0–2200 cm1), the nongray gaseous absorption coefficients due to H2O, CO2, O3, CH4, N2O, and CFC based on the correlated k-distribution method developed by Fu and Liou [15] are incorporated in the present Monte Carlo model. The infrared spectrum is divided into 12 bands 0–280, 280–400, 400–540, 540–670, 670–800, 800–980, 980–1100, 1100–1250, 1250–1400, 1400–1700, 1700–1900, and 1900–2200 cm1. The continuum absorption of H2O is included in the spectral region 280–1250 cm1. The gaseous optical depth was selected from 108 to 103. Consider a slab gas at constant temperature and constant pressure, the gaseous emissivity is defined by eu ¼ 1  T fn ðtÞ, T fn ðtÞ

is the slab transmittance. 8 6 P > > > 1  expð ei ti Þ; > > > i¼0 > > > < 6 P i e ¼ 1  expð ei t Þ; > i¼0 > > > > 6 P > i > > > 1  expð ei ðlog10 tÞ Þ; : i¼0

where

(22) The cubic gas cell emissivity can be parameterized as follows: to104 ; 104 pto1:0;

(23)

tX1:0;

where coefficients ei are obtained by numerical fitting to the emissivities computed from the Monte-Carlo model. For a homogeneous and isothermal cubic cell defined by the phase function, the single-scattering albedo, and the extinction coefficient, M photons are emitted from one of the six sides. The Monte Carlo method is then used to keep track of the individual photon, which undergoes scattering and absorption in the cell. The path of this photon will come to an end when it is either absorbed P or leaves the cubic cell. After tracking M photons, the absorptivity for the cubic cell is computed by A ¼ 1  6i¼1 M i =M, where Mi denotes the number of photons that escape from the cell surfaces. Under thermodynamic equilibrium, emissivity, e ¼ A. The fitting for cubic gaseous cell emissivity as a function of optical depth is shown in Fig. 1(a). Because the scattering effects of gases can be neglected in the thermal infrared region, the only factor affecting the gaseous absorption is its optical depth. Thus, we may parameterize the cubic gaseous emissivity as a function of optical depth for the entire thermal infrared region. This figure shows that the emissivities computed from parameterization fit closely to the ‘‘exact’’ values, especially for optical depths larger than 104. Here, the ‘‘exact’’ values refer to those computed from the Monte Carlo method. We shall use ice clouds as an example in the cloud emissivity parameterization. The single-scattering properties of ice clouds can be parameterized in terms of ice water content, IWC ,and the mean effective ice crystal size, Dge, as follows [18]: .  (24) bext ¼ IWCða0 þ a1 Dge þ a2 D2ge Þ, ba ¼

IWC ðb0 þ b1 Dge þ b2 D2ge þ b3 D3ge Þ, Dge

(25)

g ¼ c0 þ c1 Dge þ c2 D2ge þ c3 D3ge ,

(26)

$ ¼ 1  ba =bext ,

(27)

where the coefficients ai, bi, and ci are functions of wavelength. The radiative properties of ice clouds in the infrared are dominated by absorption although the scattering effects can not be neglected [19]. For this reason, the ice cloud emissivity may be directly related to the

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

171

Emissivity

1.0 0.8 0.6

Emissivity

0.4 100

102

0.2 10-3

10-2

10-1

100

10-7

10-6 Vertical Optical Depth

10-5

10-4

10-6 10-8 10-8

(a)

1.0

1.0

2200 -1900 cm-1

0.8 Emissivity

Emissivity

0.8 0.6 0.4 0.2 0

1.0

0.4

1.0

1700 -1400 cm-1 Emissivity

0.4 0.2

(b)

0.6

0

10-2 10 -1 100 101 Absorption Optical Depth

0.6

0

1900 -1700 cm-1

0.2

0.8 Emissivity

103

0.4

0.0 10-4 Emissivity

101

10-2 10 -1 100 101 Absorption Optical Depth 1400 -1250 cm-1

0.8 0.6 0.4 0.2

10 -2 10-1 100 101 Absorption Optical Depth

0

10-2 10-1 100 101 Absorption Optical Depth

Fig. 1. (a) The cubic gas cell emissivity as a function of optical depth for tX1.0 (upper panel), 104pto1.0 (middle panel), and to1.04 (bottom panel). (b) The cubic cloud cell emissivity as a function of absorption optical depth for the spectral intervals 2200–1900 cm1 (upper left panel), 1900–1700 cm1 (upper right panel), 1700–1400 cm1 (lower left panel), and 1400–1250 cm1 (lower right panel).

absorption part of optical depth (ta ¼ ba dz) in the form 8 > < 1  expð1:66ta Þ; ta o0:5; 6 P e¼ 1  expð d i tia Þ; ta X0:5; > : i¼0

(28)

ARTICLE IN PRESS 172

Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

for PP and homogeneous ice clouds and ! 6 X i e ¼ 1  exp  f i ta

(29)

i¼0

for cubic ice clouds, where the coefficients d i and f i are obtained by numerically fitting of the emissivities computed from the Monte Carlo model for six ice crystal size distributions [20,21]. The IWC ranges from 0.0002 to 0.2 g m3. Fig. 1(b) shows the fitting for the cubic ice-cloud cell emissivity as a function of absorption optical depth in the spectral intervals 2200–1900 cm1 (upper left panel), 1900–1700 cm1 (upper right panel), 1700–1400 cm1 (lower left panel), and 1400–1250 cm1 (lower right panel). For the same absorption optical depth, the ice cloud emissivities differ for different spectral intervals, due primarily to the scattering effects of ice particles, which are spectral dependence. 2.3. Monte Carlo statistical error The uncertainty in flux quantities cannot be computed directly from the square root of the number of photons as is usually done when nonweighted Monte Carlo techniques are used [2]. The ability to calculate uncertainties is crucial when the accuracy of model results is dependent on the number of photons used in the Monte Carlo technique. We have used a simple technique [22] to compute the upper limit on the uncertainty in flux quantities for each cell as follows: PP m 1 sX ði; j; kÞ2 ¼ Mði;j;kÞðMði;j;kÞ1Þ ðxn  X ði; j; kÞÞ2 m n (30) 2 1 ¼ Mði;j;kÞ1 ðX ði; j; kÞ  ðX ði; j; kÞ2 Þ; where Mði; j; kÞ is the number of the scattering event that contributes to the output quantity in the cell surface ði; j; kÞ. The sums over m and n were done only for those photons that contribute to the flux quantity in the cell surface ði; j; kÞ, while X ði; j; kÞ is calculated by using Eqs. (16)–(21). The simplest and probably the most pragmatic approach is to judge the accuracy of model output by monitoring its convergence and stability as a function of the number of photons used in the simulation. Fig. 2 30 12,000K photon count 24,000K photon count 36,000K photon count 48,000K photon count

20

% of flux difference

10

0

-10

-20

-30

-40 0

10

20

30 40 Height (km)

50

60

70

Fig. 2. Relative flux differences at several heights for different photon counts. The solid lines display the upward flux relative difference to the 96,000 K photon count at each height, while the dash lines show the downward flux relative difference.

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

173

shows that the upward flux difference is within 0.1% at all heights regardless of the photon number, while the downward flux difference varies significantly with the photon number and height. The downward flux difference is within 1% below 30 km, but can be as large as 40% at TOA because of the lower air density. Relative accurate fluxes can be obtained by using 50,000 K photons up to 50 km. 2.4. Parallelization of the Monte Carlo thermal emission model In order to achieve convergence with a minimal statistical fluctuation, a large number of photons are required. The present 3D Monte Carlo model, which loops over the grid cell, is extremely time-consuming. Thus, we developed a vectorized program to speed up the calculation. The resulting vectorized program was run on a super-computer at NCAR. Ray tracings of the multiple photons emitted from the surfaces of each grid cell are performed for the 12 thermal infrared bands each with a different number of weights determined by the correlated-k distribution method. Considering all the weighting points, a total of 67 monochromatic bands is required to obtain broadband flux quantities. We may divide a portion of the bands among the computer processors. Each processor would then determine the number of rays originated from a given cell in their domain to be traced. The rays are then traced following the conventional technique. This method is effective because it involves less communication among photons and only requires to past the band number to other processors. We used the message passing interface to implement the parallel version of the model because of the ease with which it could be incorporated into the serial model framework. The master thread packs a band number to the slaves. After sending the slaves their first ‘‘band’’, the master waits for one to signal that it is finished, at which point the master sends the next band number off. This continues until the entire band has been parceled out. At the same time, the master signals to the slaves that they are done. The flux quantities are store in temporal files by the slave processors. The master processor sums them up to obtain the broadband infrared flux quantities. 3. Computational results Since the present 3D thermal infrared Monte Carlo (MC) method is new and no previous results are available for cross check, we compare our calculations to those obtained from the PP model under the same input condition. The well developed Fu–Liou model [8,9,15,16] for broadband radiative transfer calculations in plane-parallel atmospheres has been used to generate pertinent results for comparison. This model has been validated via comparison to an ‘‘exact’’ line-by-line method for a variety of clear atmospheres using the same input. Tables 1 and 2 list the upward flux at TOA and the downward flux at the surface for the MC and Fu–Liou PP models. We have carried out two types of MC calculations. One involves the PP condition, that

Table 1 Comparison of the Fu–Liou (FL) and 1D and 3D Monte Carlo (MC) radiative transfer schemes for broadband infrared fluxes (W m2) at the top of the atmosphere (TOA) and the surface (SFC) using different atmospheric profiles under the clear sky condition. 30 million photons are used in the MC simulation Atmospheric conditions

F " ðTOAÞ FL

MLS MLW SAS SAW TRO USA

284.2 233.1 265.7 200.7 292.1 262.9

F # ðSFCÞ MC

FL

1D

3D

283.3 232.5 264.6 200.1 292.2 261.9

284.9 233.4 266.1 201.1 292.8 263.1

348.6 221.2 299.1 170.3 396.1 284.6

MC 1D

3D

346.2 220.9 296.5 171.8 393.5 282.1

346.7 220.3 296.4 171.6 393.5 281.6

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

174

Table 2 Same as in Table 1, except a high cloud ðt ¼ 4; De ¼ 42mmÞ is placed between 9 and 11 km Atmospheric conditions

F " ðTOAÞ FL

MLS MLW SAS SAW TRO USA

F # ðSFCÞ MC

173.1 138.4 156.6 128.4 175.4 147.0

FL

1D

3D

170.6 135.7 153.8 125.8 173.1 144.2

170.3 135.5 153.7 125.9 172.4 144.1

50 MC MC Fu-Liou

30 20

40

(a)

-8 -6 -4 -2 Cooling Rate (K/day)

0 -0.2 0 0.2 0.4 0.6 Difference (K/day)

MC

20

40 Height (km)

Height (km)

30

-6 -2 Cooling Rate (K/day)

(b)

20 10

(c)

366.7 247.6 317.3 199.6 407.3 306.5

0 -0.6 -0.4 -0.2 0 0.2 Difference (K/day)

50

30

0 -12

367.0 248.8 318.0 200.3 408.1 307.6

MC

0 -10

50 40

3D

10

10 0 -10

1D

50

Height (km)

Height (km)

40

368.9 248.9 320.2 198.9 410.5 309.9

MC

MC

30 20 10

-8 -4 Cooling Rate (K/day)

0

-0.6 -0.4 -0.2 0 0.2 Difference (K/day)

0 -12 (d)

-8 -4 Cooling Rate (K/day)

0 -0.6 -0.4 -0.2 0 0.2 Difference (K/day)

Fig. 3. The infrared broadband cooling rate profiles computed from the 3D Monte Carlo method (solid line) and differences (dot line) from the Fu–Liou model in clear sky. (a) The subarctic winter atmosphere (SAW), (b) the tropical atmosphere (TRO), (c) the mid latitude summer atmosphere (MLS), and (d) the 1976 US standard atmosphere (USA).

is, variation of the atmospheric parameters is only in the vertical (1D). The other uses the cubic cell as the building block (3D) but assuming that all the cubes have the same optical properties horizontally. For 3D MC, the boundary condition used was cyclic in the horizontal plane, i.e., when the photon leaves the domain from one side; it will enter into the domain from the opposite side. Six atmospheric profiles, including the midlatitude summer (MLS), midlatitude winter (MLW), subarctic summer (SAS), subarctic winter (SAW), tropical atmosphere (TRO), and the 1976 US Standard Atmosphere (USA) were used in flux and heating rate

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

175

calculations. The CO2, CH4, and N2O concentrations were set as 370, 1.6, and 0.28 ppmv, respectively. In the calculations, the atmosphere was equally divided into a number of homogeneous layers with a vertical resolution of 0.5 km and the surface emissivity was set to 1. For the cloudy case, a high cloud with an optical depth of 4 and a mean effective ice crystal size of 42 mm was placed between 9 and 11 km. Results computed from the 1D MC model are in excellent agreement with those from the Fu–Liou PP model. For the upward flux at TOA, differences are generally less than 1 W m2, while for the downward flux at the surface, differences are within about 2 W m2. The flux differences between 1D and 3D MC calculations defined above are generally within 0.5 W m2. The slight differences are produced by parameterization of the emissivities for gases [Eqs. (22) and (23)] and clouds [Eqs. (28) and (29)] in 1D and 3D cases. Fig. 3 shows the infrared broadband heating rate profiles for the 3D MC and the heating rate difference between 3D MC and Fu–Liou models in clear sky. The differences in broadband upward, downward, and net fluxes between the two (not shown) are less than 1.0 W m2 for all layers in clear sky. The absolute differences in heating rates are less than 0.1 K day1 for all layers below 30 km. But above about 30 km larger differences are shown due to the smaller emissivity for gases and the consequence of a smaller number of emitted photons, as compared with lower layers in the Monte Carlo model. For cloud sky (Fig. 4), differences between the two models for fluxes and cooling rates are on the same order as in clear sky conditions. However, larger differences of about 0.4 K day1 occur in the cloud layer. Comparisons of the MC results for broadband thermal cooling rates in 3D clouds to those computed from the delta-diffusion approximation (DDA) for 3D radiative transfer [9] and the independent pixel-by-pixel approximation (IPA) are subsequently carried out to understand the relative merits of these approaches. For these comparisons, we selected the optical depth and mean effective ice crystal size over an area of 30 km by

50

50 MC

MC MC Fu-Liou

30 20

40 Height (km)

Height (km)

40

(a)

20 10

10 0 -30

30

-20 -10 0 10 Cooling Rate (K/day)

20 -0.2 0 0.2 0.4 0.6 Difference (K/day)

0 -30 (b)

50 MC

MC

Height (km)

Height (km)

40

30 20 10

(c)

20 -0.6 -0.4 -0.2 0 0.2 0.4 Difference (K/day)

50

40

0 -30

-20 -10 0 10 Cooling Rate (K/day)

30 20 10

-20 -10 0 10 Cooling Rate (K/day)

20 -0.6 -0.4 -0.2 0 0.2 0.4 Difference (K/day)

0 -20 (d)

-10 0 10 Cooling Rate (K/day)

20

-0.2 0 0.2 0.4 Difference (K/day)

Fig. 4. Same as in Fig. 3, except a high cloud (t ¼ 4,De ¼ 42 mm) is placed between 9 and 11 km.

ARTICLE IN PRESS 176

Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

Fig. 5. Comparison results for the domain averaged broadband infrared cooling rates. (a) Results computed from the 3D Monte Carlo method, (b) relative difference (%) between the 3D delta-diffusion model and the Monte Carlo method, and (c) relative difference (%) between the independent pixel-by-pixel approximation and the Monte Carlo method. The 3D results are averaged over the height and presented in the horizontal plane (upper panel) and averaged over the x direction and shown in the yz-plane (lower panel).

20 km determined from the NOAA advanced very high resolution radiometer (AVHRR) presented in Gu and Liou [8] and Chen et al. [9]. Fig. 5(a) displays the 3D MC results presented in the horizontal plane (averaged over the height), and in the yz-plane (averaged over the x-direction). Figs. 5(b) and (c) depict the relative differences for 3D DDA and IPA with respect to the 3D MC method, respectively. In the horizontal plane (xyplane) the relative differences range from 10% to 10% for DDA and IPA, respectively. For smaller IWCs, the 3D DDA produces smaller cooling rates (negative relative differences), whereas for larger IWCs, the opposite is true. The 3D DDA approach appears to have more horizontal flux transfers among cloud columns, as compared to the 3D MC simulation. The situation is reverse for the IPA results. For smaller IWCs, IPA generates relative larger cooling rates (positive relative differences), whereas relative smaller cooling rates are produced for larger IWCs (negative relative differences) due to the fact that radiative fluxes are computed on a pixel-by-pixel basis without horizontal flux exchange between pixels. The pronounced heating at the cloud bottom and the strong cooling at the cloud top are produced by the absorption of cloud particles and the subsequent emission of these particles in the corresponding infrared region associated with cloud temperature. The net horizontal radiative transfer is very important to the flux and cooling rate distributions for inhomogeneous clouds. For the vertical plane (yz-plane) the relative differences range from 10% to 2%. The relative differences between 3D DDA and 3D MC are very small, less than 4%. The maximum difference occurs in the middle of the cloud layer. However, IPA underestimates the heating rate at the cloud bottom by about 10%, and overestimates the cooling rate at the cloud top by about 2%. Finally, it is noted that the 3D MC method requires much more computer time (by a factor of 100) as compared to DDA and IPA. 4. Summary We have developed a 3D broadband thermal radiative transfer scheme by integrating a 3D thermal infrared Monte Carlo photon transport algorithm into the Fu–Liou spectral radiation program. In this model the gaseous absorption and cloud emissivity are parameterized by using a single cubic element as the building block, which can be applied to radiative transfer in complex cloud geometry. The validity of the present Monte Carlo approach was examined by comparing its 1D result to the plane-parallel Fu–Liou counterpart under a

ARTICLE IN PRESS Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

177

variety of atmospheric conditions. We found that the results computed from the Monte-Carlo model are in excellent agreement with those from the Fu–Liou plane-parallel model. Differences are generally less than 1 W m2 for the upward flux at TOA, while for the downward flux at the surface, differences are within about 2 W m2. The flux differences between 1D and 3D MC calculations are generally within 0.5 W m2. The absolute differences in heating rates are less than 0.1 K day1 for all layers below 30 km under clear sky. The largest differences for cloud layers are small than 0.4 K day1. We further examined the present Monte Carlo approach for thermal infrared radiative transfer in 3D clouds by comparing to the results computed from the delta-diffusion approximation for 3D radiative transfer and the independent pixel-by-pixel approximation, which applies the plane-parallel radiation calculation to each column of the cloud field. Relative differences are found to be within about 10% for the domain-averaged heating rates. The comparison results also illustrate that the net horizontal radiative transfer is important in the calculation of flux and cooling-rate distributions for inhomogeneous clouds. Finally, it should be noted that although the present 3D Monte Carlo method for thermal infrared radiative transfer requires substantial computer time to achieve reliable accuracy for flux calculations, it can be used as a benchmark for the development and validation of other analytical 3D thermal radiative emission approaches. Acknowledgements This work has been supported by DOE Grant DE-FG02-04ER6333724, NSF Grant ATM-0331550, and NASA Grant NNG04GG91G.

References [1] Plass GN, Kattawar GW. Monte Carlo calculations of light scattering from clouds. Appl Opt 1968;7:415–9. [2] Cahalan R, Ridgway W, Wiscombe WJ, Gollmer S, Harshvardhan. Independent pixel and Monte Carlo estimates of stratocumulus albedo. J Atmos Sci 1994;51:3776–90. [3] O’Hirok W, Gautier C. A three-dimensional radiative transfer model to investigate the solar radiation within a cloudy atmosphere. Part I: spatial effects. J Atmos Sci 1998;55:2162–79. [4] Barker HW, Morcrette JJ, Alexander GD. Broadband solar fluxes and heating rates for atmospheres with 3D broken clouds. Q J R Meteorol Soc 1998;124:1245–71. [5] Fu Q, Cribb MC, Barker HW, Krueger SK, Grossman A. Cloud geometry effects on atmospheric solar absorption. J Atmos Sci 2000;57:1156–68. [6] Liou KN, Ou SC. Infrared radiative transfer in finite cloud layers. J Atmos Sci 1979;36:1985–96. [7] Evans KF. The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer. J Atmos Sci 1998;55:429–46. [8] Gu Y, Liou KN. Radiation parameterization for three-dimensional inhomogeneous cirrus clouds: Application to climate models. J Clim 2001;14:2443–57. [9] Chen Y, Liou KN, Gu Y. An efficient diffusion approximation for 3D radiative transfer parameterization: application to cloudy atmospheres. JQSRT 2005;92:189–200. [10] Gentile N A. Implicit Monte Carlo diffusion—an acceleration method for Monte Carlo time-dependent radiative transfer simulations. J Comput Phys 2001;172:543–71. [11] Takara EE, Ellingson RG. Broken cloud field longwave-scattering effects. J Atmos Sci 2000;57:1298–310. [12] Liu Q, Simmer C, Ruprecht E. Three-dimensional radiative transfer effects of clouds in the microwave spectral range. J Geophys Res 1996;101(D2):4289–98. [13] Roberti L, Kummerow C. Monte Carlo calculations of polarized microwave radiation emerging from cloud structures. J Geophys Res 1999;104(D2):2093–104. [14] Davis C, Emde C, Harwood R. A 3D polarized reversed Monte Carlo radiative transfer model for mm and sub-mm passive remote sensing in cloudy atmospheres. IEEE TGARS, MicroRad04 Special Issue, 2004. [15] Fu Q, Liou KN. On the correlated k-distribution method for radiative transfer in nonhomogeneous atmospheres. J Atmos Sci 1992;49:2139–56. [16] Fu Q, Liou KN. Parameterization of the radiative properties of cirrus clouds. J Atmos Sci 1993;50:2008–25. [17] Liou KN. An Introduction to Atmospheric Radiation, 2nd ed. San Diego: Academic Press; 2002. [18] Fu Q, Yang P, Sun WB. An accurate parameterization of the infrared radiative properties of cirrus clouds for climate models. J Clim 1998;11:2223–37. [19] Fu Q, Liou KN, Cribb MC, Charlock TP, Grossman A. Multiple scattering parameterization in thermal infrared radiative transfer. J Atmos Sci 1997;54:2799–812.

ARTICLE IN PRESS 178

Y. Chen, K.N. Liou / Journal of Quantitative Spectroscopy & Radiative Transfer 101 (2006) 166–178

[20] Heymsfield AJ, Platt CMR. A parameterization of the particle size spectrum of ice clouds in terms of the ambient temperature and ice water content. J Atmos Sci 1984;41:846–55. [21] Takano Y, Liou KN. Solar radiative transfer in cirrus clouds. Part I: single scattering and optical properties of hexagonal ice crystals. J Atmos Sci 1989;46:3–19. [22] Gordon KD, Misselt KA, Witt AN, Clayton GC. The dirty model. I. Monte Carlo radiative transfer through dust. Astro J 2001;551:269–76.