## Abstract

A new calculation scheme is proposed for the explicitly discretized solution of the three-dimensional (3D) radiation transfer equation (RTE) for inhomogeneous atmospheres. To separate the independent variables involved in the 3D RTE approach, the spherical harmonic series expansion was used to discretize the terms, depending on the direction of the radiance, and the finite-volume method was applied to discretize the terms, depending on the spatial coordinates. A bidirectional upwind difference scheme, which is a specialized scheme for the discretization of the partial differential terms in the spherical harmonic-transformed RTE, was developed to make the equation determinate. The 3D RTE can be formulated as a simultaneous linear equation, which is expressed in the form of a vector–matrix equation with a sparse matrix. The successive overrelaxation method was applied to solve this equation. Radiative transfer calculations of the solar radiation in two-dimensional cloud models have shown that this method can properly simulate the radiation field in inhomogeneous clouds. A comparison of the results obtained using this method with those using the Monte Carlo method shows reasonable agreement for the upward flux, the total downward flux, and the intensities of radiance.

## 1. Introduction

Analytic methods commonly employ discretization schemes to separate the independent variables of a 3D RTE. For example, Stephens (1988a) and Gabriel et al. (1993) employed the Fourier series expansion to discretize the terms that depend on the horizontal coordinates, as well as the Gauss quadrature that was used to represent the direction of radiance. Ishida and Asano (2007) used a similar method but used a scaling function expansion for the discretization. This scheme is effective in expressing the fine horizontal structures of the radiation fields with smaller terms. In these methods, the 3D RTE consequently results in a vector–matrix differential equation that is analogous to the 1D RTE in a plane-parallel atmosphere. The only independent variable of the equation is the vertical coordinate. Thus, a method used for the 1D RTE solution, for example, the doubling–adding method (Stephens 1988a; Ishida and Asano 2007) or the eigenvalue method (Kobayashi 1991), can be applied to solve the boundary problem of the 3D RTE. On the other hand, Evans (1993) proposed the Spherical Harmonic Spatial Grid (SHSG) method that divides the space into many infinitesimally small volumes to discretize the terms that depend on the horizontal coordinates and also on the vertical coordinate. The SHSG method also uses a spherical harmonics series to separate the independent angular variables for the direction of radiance; that is, the zenith and azimuth angles. Later, Evans (1998) provided the Spherical Harmonics Discrete Ordinate Method (SHDOM), which combines the spherical harmonics functions and the discrete ordinate method to express the direction of radiance. The SHDOM is versatile for any radiation field, but the transform between the spherical harmonics functions and the discrete ordinate method is a complicated procedure.

This paper provides a new calculation scheme for an explicitly discretized solution to the 3D RTE. This method is an improvement on the SHSG method, but it is more efficient and suitable for fine and anisotropic radiation fields, such as solar radiation in optically thin clouds. This method uses the spherical harmonic expansion to discretize the terms that depend on the direction of radiance, as it facilitates calculations of the radiance in an arbitrary direction. Our method also employs the finite-volume method to discretize the terms that depend on the spatial coordinates. We have modified the upwind difference scheme used for the expression of the spatial differentiation to make the transformed RTE determinate, because the general upwind difference scheme is not appropriate for the transformed RTE, as discussed below. In section 2, we describe the algorithm of the solution method in detail. In section 3, we present calculation examples for simple cloud models and validate our method by comparing it to data obtained using the Monte Carlo method. The dependence of the accuracy of the discretization on the spatial and angular resolution is discussed in section 4, and in section 5, we present a summary and our conclusions.

## 2. Description of the algorithm

### a. Discretization of the angle-dependent terms

In this section, we formulate the discretized 3D RTE in a vector–matrix form. A general 3D RTE for monochromatic radiance is expressed as follows (Ishida and Asano 2007):

where the definitions of the variables are listed in Table 1. The left-hand side of Eq. (1) is known as the “streaming term” (Marshak and Davis 2005). Such terms depend on the azimuth and zenith angles, as N, p, and Q can be discretized in the forms of a spherical harmonics series. For example, the radiance, N(r; μ, ϕ), is expressed in the form

where Nc/sl,m is the transformed intensity with spherical harmonics functions, Yc/sl,m(μ, ϕ), that are defined by

where Pml(μ) is the associated Legendre function of degree = l and order = m, and L and M are the maximum number of degrees and orders for the expansion, respectively. Here, M must be equal to or smaller than L.

Table 1.

Variables used in the three-dimensional radiative transfer equation.

The spherical harmonic transform of degree, l, and order, m, for the streaming term is determined from the transformed radiances using the spherical harmonic functions of the other degree and order numbers. For example, the transform of sinθ cosϕ(∂N/∂x) with Yc/sl,m(μ, ϕ) results in

where

and the al,m are zero for values of l′ and m′ other than l ± 1 and m ± 1. When l′ or m′ is above the maximum number (L or M), then the al,m are also set to zero. Thus, the streaming term is given by

where

The phase function, p(r; μ, ϕ; μ′, ϕ′), can be expanded in a series of Legendre polynomials, Pl(cosΘ), in the form

where χl is the l-degree coefficient of the Legendre expansion of the phase function and Θ is the scattering angle. Using the addition theorem, the spherical harmonic transform of the scattering integral term results in the following simple form (Evans 1993):

After the spherical harmonic expansion, Eq. (1) is finally expressed in the form

Furthermore, Eq. (11) can be written as a vector–matrix equation, that is,

where the notation of the matrices follows that of Marshak and Davis (2005). The elements of matrices 𝗠x, 𝗠y, and 𝗠z are al,m, bl,mc/s, and cl,m, respectively, and the elements of matrix 𝗖 are α + sχl/(2l + 1), and vector S consists of the transformed source terms, which are the 2nd and 3rd terms of the right-hand side of Eq. (11).

### b. Discretization of the spatial coordinates

The terms that depend on the horizontal coordinates are discretized using the finite-volume method, which has been applied to obtain solutions of differential equations in many fields, such as heat transfer (e.g., Chai et al. 2004). Discretization using the finite-volume method satisfies the conservation law of physical quantities (e.g., energy) both locally and fully. The finite-volume method requires the use of a difference scheme to represent the partial differential terms. There are several types of difference schemes, depending on the number of grid points used and the precision required. For example, the central difference scheme, which was used in the SHSG method, is applicable to an equation that involves small advection terms (which correspond to the streaming term in the RTE), or the large diffusion terms. However, since the 3D RTE does not involve any diffusion terms, the use of the central difference scheme for solving the 3D RTE sometimes causes unrealistic numerical oscillations in the radiation field, or no convergence of the solution occurs. In particular, the discretized 3D RTE for a cloud layer that contains vacant or no absorption areas results in an indeterminate equation. The central difference scheme is also not appropriate to a boundary value problem, such as the 3D RTE, because the treatment of boundary conditions becomes complicated. On the other hand, the upwind difference scheme is suitable for discretization of an equation that involves large advection terms. When a scalar variable moves as a flow, the upwind difference scheme uses only the quantities on the upward side of the flow to represent the difference of a scalar variable. However, the “velocity” of a radiance flow at a given position cannot be determined uniquely, because the radiance itself involves direction and can come from any direction. Therefore, we had to modify the upwind difference scheme for adapting to such bidirectional flows.

In this section, we will explain the finite-volume method with the bidirectional upwind difference scheme for a 3D radiation field. A schematic of the finite-volume coordination is shown in Fig. 1. First, we consider a cubic cell whose center is located at a grid point P in Figs. 1a,b. The length of the cell along the x-, y-, and z-axis coordinates is denoted as Δx, Δy, and Δz, respectively. It is necessary to determine the spherical harmonic-transformed radiance, Nc/sl,m, on each boundary of the cell (e, w, b, t in Fig. 1b). The Nc/sl,m on a boundary can be divided into two components of the positive and negative directions along the coordinates by applying the hemispherical integral. For example, on the boundary that is parallel to the zy plane (corresponding to the plane including position e in Fig. 1b), Nc/sl,m can be written as

where J+e;c/s,l,m and Je;c/s,l,m are the components of the spherical harmonic-transformed radiances for the positive direction (rightward direction in Fig. 1b) and the negative direction (leftward direction) along the x coordinate at position e. The terms, J+e;c/s,l,m and Je;c/s,l,m, can be expressed by

and

where Rc/s;l,m;l′,mrgt and Rc/s;l,m;l′,mlft are the coefficients needed for obtaining the hemispherical integral with the spherical harmonic functions. The detailed definitions of these terms are given in appendix A.

Fig. 1.

A schematic drawing of the finite-volume method: (a) the symbols of the neighborhood grid points of position P and the virtual cell that includes position P, and (b) the symbols of the boundaries of the cell (in two dimensions) and the range of the hemispherical integral on a boundary. Case (i) is a general cell in the domain. The terms Je+ and Je denote the forward and backward directions of radiance along the x coordinate on the boundary, e. Case (ii) is a cell for a vertical boundary of the domain and case (iii) is for a horizontal boundary.

Fig. 1.

A schematic drawing of the finite-volume method: (a) the symbols of the neighborhood grid points of position P and the virtual cell that includes position P, and (b) the symbols of the boundaries of the cell (in two dimensions) and the range of the hemispherical integral on a boundary. Case (i) is a general cell in the domain. The terms Je+ and Je denote the forward and backward directions of radiance along the x coordinate on the boundary, e. Case (ii) is a cell for a vertical boundary of the domain and case (iii) is for a horizontal boundary.

The upwind difference scheme assumes that a physical quantity at a boundary is the same as that of the upwind grid point. For example, J+e;c/s,l,m is determined from the (transformed) radiances at grid point P (i.e., the upwind position of the positive direction along the x coordinate at position e) in the form

where rP is the position vector at grid point P. On the other hand, Je;c/s,l,m is determined from the radiance at the grid point E (i.e., the upwind position of the negative direction along the x coordinate at the position e) as

where rE is the position vector at the grid point E (determining J+e;c/s,l,m from the mean value of the grid points P and E corresponding to the center difference scheme). Equations (16a) and (16b) can be written in the vector–matrix form

and

In a similar way, the components of the spherical harmonic-transformed radiances for the positive and negative direction along the x coordinate at position w, J+w and Jw, are given by

and

where rW is the position vector at the grid point W. The corresponding quantities along the y coordinate at position n (i.e., front of the cell) and position s (i.e., rear of the cell), J+n, Jn, J+s, and Js, are given by

where rN and rS are the position vectors at grid points N and S, respectively. Similarly, those along the z coordinate at position t (i.e., top of the cell) and position b (i.e., bottom of the cell) in Fig. 1b, J+t, Jt, J+b, and Jb, are given by

where rT and rB are the position vectors at grid points T and B, respectively. The elements of the matrices, Rc/s;l,m;l′,mfrn, Rc/s;l,m;l′,mbck, Rc/s;l,m;l′,mup, and Rc/s;l,m;l′,mdwn, are also given in appendix A.

Next, we take into account the conservation of a physical quantity (corresponding here to Nc/sl,m) in the cell, which is considered as a finite volume. We assume that the cell is small enough so that a physical quantity can be considered to be constant on each of the cell boundaries. From this assumption, the net outflow of a physical quantity through a boundary can be easily determined from the volume integral of the left-hand side of Eq. (12), and this volume integral can be transformed to the surface integral using Gauss theorem. For example, the net outflow through plane e is

The source of a physical quantity in the cell can also be estimated by the volume integral of the right-hand side of Eq. (12). Assuming that the right-hand side of Eq. (12) is uniform in the cell, the source is

Since the source of a physical quantity must be equal to its net outflow to satisfy the conservation law of energy within the cell, Eq. (12) can be discretized in the form

In the same manner, this discretization scheme can be also applied to a cell at the domain boundary. An example is cell (ii) in Fig. 1b, where the net flow at the cell boundary is Pb and results in

where

are the components of the spherical harmonic-transformed radiances coming in from the outer environment of the domain through the boundary. Use of the radiance at the opposite side of the domain for the incoming radiance through the boundary corresponds to the periodic boundary condition. The horizontal boundary condition can be defined in a similar way. In particular, reflecting the surface boundary condition for the grid point Pc that is at the bottom boundary, corresponding to cell (iii) in Fig. 1b, is given by

where matrix 𝗥s represents reflection of the surface and its elements are defined in appendix B.

The direct incident radiance is calculated from the integral of the extinction coefficients along the traveling path of the radiance. For this integral, we can apply a cubic spline interpolation of the extinction coefficients from the values at the grid points. Thus, it should be noted that the direct component of the incident solar radiation could be obtained without solving the RTE.

### c. The solution method for the vector–matrix equation

Since the set of Eq. (25) for all positions results in a simultaneous linear equation, this can be expressed in a simple vector–matrix form:

where N is the vector that contains all the terms of the spherical harmonic series at all grid points, and e is the corresponding vector that consists of the discretized source terms. Matrix 𝗔 is very sparse because of the characteristics of the spherical harmonic transform, as in Eq. (10). The physical properties, such as the extinction coefficient and the phase function, can be independently defined for each grid point. There are several methods that can be used to solve Eq. (29). A most primitive method is to obtain the inversion matrix of 𝗔, which requires a large effort, but can explicitly represent the relationship between a cloud field and the corresponding radiation field. Therefore, this method may be helpful in investigating the effect of cloud inhomogeneity on radiation. Another method is the conjugate gradient method (Evans 1993), which estimates an approximation value using iterative calculations, and is versatile for solving a simultaneous linear equation system. In this study, we employed the successive overrelaxation (SOR) method, which is a relaxation method that is often used to solve a differential equation for a steady-state flow. For example, Gu and Liou (2001) applied the SOR method to solve the RTE. This method is appropriate for solving a difference equation that consists of values at a grid point and its neighbors, such as in Eq. (25), and enables us to obtain the radiances at all of the grid points in a considered domain simultaneously. An advantage of the SOR method is that it only requires a single large array for storing the values of the discretized radiances. The SOR method iterates the following calculation,

until the values converge, where matrices 𝗗, 𝗟, and 𝗨 represent the diagonal, strictly lower-triangular, and strictly upper-triangular parts of 𝗔, respectively, and r is referred to as the relaxation parameter. For the initial radiance values, we selected the source function values obtained from the first single scattering of the incident solar radiation. In general, a larger relaxation parameter requires fewer iterations. We set the relaxation parameter to a value of 0.5, which lets all the difference equations have convergent solutions, because our calculations on the models described below with a relaxation parameter greater than 0.5 sometimes failed to obtain a convergent solution. This may imply poor convergence properties of the discretized 3D RTE, especially for cases with a vacant cloud field and/or conservative scattering. The criterion for convergence was set to satisfy |N(r)n+1N(r)n|/|N(r)n| < 10−4. One step of the iteration of the SOR method takes one matrix–vector multiplication, two vector norm calculations, two vector–scalar multiplications, and two vector additions. When L is odd and M is even, the required number of floating point operations per one iteration is 2n2 + {2(L + 1)(M + 2) + 2(L + 3) + L + 3 + (L + 1)(M/2 + 1) + 3} × n, where n = (L + 1) × (M + 1) × nx × ny × nz and nx × ny × nz is the number of spatial grid points. After N in Eq. (29) is determined, the intensities of scattered radiance at a grid point can be determined by converting from the spherical harmonic series as in Eq. (2).

## 3. Examples and validation

### a. The cloud model

In this section, we will validate the performance of our solution by showing calculations of radiation fields using hypothetical cloud models, and we will show a comparison of the data with that obtained using the Monte Carlo method. The influence of the angular and spatial resolutions on the calculation accuracy was also investigated. For simplicity, we considered a rectangular-shaped inhomogeneous cloud field. The origin of the Cartesian coordinates was set in the bottom-left corner of the cross section of the domain. The domain was assumed to be uniform and to extend in an infinite direction along the y coordinate; that is, the cloud model is two dimensional along the x and z coordinates. The geometrical thickness of the domain was set to 800 m, and the length along the x axis was 3.2 km. We only considered the transfer of the monochromatic solar radiation at a wavelength of 500 nm and neglected any molecular scattering and absorption. The solar incident angle was set to 53°, with a cosine of 0.6. The single scattering albedo of the cloud particles was assumed to be unity, which implied the conservative scattering case. We assumed a uniform cloud-particle size distribution, but with different concentrations throughout the entire domain. The size distribution obeyed a lognormal function, with an effective radius of 4 μm and a standard deviation of 1.5. We assumed a water cloud, and the scattering phase function of the cloud particles was calculated using Mie’s theory. The phase function was truncated using the delta-M method (Wiscombe 1977) to deal with any strong forward scattering by the cloud particles. The truncation number for the delta-M method had to be equal to, or less than, L in Eq. (2). In this study, we selected L to represent the truncation number. A smaller truncation number can lead to larger errors, especially in radiances at optically thin points. The direction of the incident solar radiation was set to be parallel to the xz plane. For the boundary condition, we assumed that there was no inflow of diffusive radiance through all of the domain boundaries, and no reflection occurred at the ground surface. Hereafter, we will use the normalized radiative fluxes and intensities relative to the incident downward solar flux.

### b. Calculated results and comparison

The spatial distributions of the extinction coefficients for two cloud models are shown in Figs. 2a,b, in which the horizontally averaged optical thickness of each cloud model was set to a value of 8. The cloud field in Fig. 2a corresponds to a section of broken stratocumulus. The cloud model is virtual and simple, but is appropriate to use in an investigation of the effect of a broken cloud layer. Figure 3 shows the spatial distribution of the upward fluxes and total (diffusive plus direct) downward fluxes computed for the cloud model shown in Fig. 2a. The maximum degree (L) and order (M) of the spherical harmonic series expansion were set to values of 15 and 8, respectively. The interval between the spatial grid points was 50 m for the x coordinate and 25 m for the z coordinate. It can be seen that the calculated radiation field does not contain any artificial oscillations. In contrast, when the center difference scheme was applied in the discretization of the differential terms (the other calculation routines were not changed), a convergent solution could not be obtained. This result indicates that the present method is capable of simulating the radiation fluxes in inhomogeneous clouds.

Fig. 2.

(a) The spatial distribution of the cloud extinction coefficients (m−1) in a inhomogeneous cloud model for the calculation example. (b) Same as in (a), but for a more complicated case. The mean optical thickness of both cloud models was set to 8.

Fig. 2.

(a) The spatial distribution of the cloud extinction coefficients (m−1) in a inhomogeneous cloud model for the calculation example. (b) Same as in (a), but for a more complicated case. The mean optical thickness of both cloud models was set to 8.

Fig. 3.

Calculated spatial distributions of (a) the upward flux and (b) the total downward flux in the domain of the cloud model shown in Fig. 2a. The cosine of the solar incident zenith angle is −0.6 (θ0 = −53°).

Fig. 3.

Calculated spatial distributions of (a) the upward flux and (b) the total downward flux in the domain of the cloud model shown in Fig. 2a. The cosine of the solar incident zenith angle is −0.6 (θ0 = −53°).

The upward flux and the downward flux at the gap between broken clouds A and B in Fig. 2a have nonzero values, and the downward fluxes beneath the optically thin positions of the sunny side of cloud B exceed a value of 1.0. These features are due to the horizontal transport of the radiation. The net horizontal diffusive flux along the x coordinate is shown in Fig. 4 to clarify the phenomenon of horizontal transportation in the cloud layer. A large amount of the radiance penetrating into the cloud escapes from the opposite lateral side of the solar incident. Reflection of the solar incident radiance on the solar-illuminated lateral side of the cloud layer was also detected. At the gap in the cloud, the horizontal net flux was small and almost zero, because of the cancellation of the rightward flux from the left cloud block (A in Fig. 2a) by the leftward flux from the right cloud block (B). The horizontal transport through the cloud sides causes a decrease in the downward flux beneath the cloud layer compared to data from calculations using the independent columnar approximation (ICA) method. This is one of the significant effects of horizontally inhomogeneous clouds. We estimated the ratio of the outgoing energy from this domain to the total incoming energy into the domain. The ratios of energy through the top, bottom, left (x = 0.0 km), and right (x = 3.2 km) boundary are about 0.24, 0.51, 0.04, and 0.21, respectively, and the ratio of the total outgoing energy is about 1.0: this proves that the conservation of radiative energy within the domain is satisfied.

Fig. 4.

The spatial distribution of the net horizontal diffusive flux using the cloud model shown in Fig. 2a. The positive value means that the flux directed along the positive direction of the x-coordinate axis is larger than that directed along the negative direction. The values are normalized to the solar incident flux.

Fig. 4.

The spatial distribution of the net horizontal diffusive flux using the cloud model shown in Fig. 2a. The positive value means that the flux directed along the positive direction of the x-coordinate axis is larger than that directed along the negative direction. The values are normalized to the solar incident flux.

Figure 5 shows the intensities of the diffusive radiance at some directions for the upward (μ = 0.95 and μ = 0.55) and downward (μ = −0.95 and μ = −0.55) in the principal plane of ϕ = 0°. This shows that our method is capable of properly simulating the nature of radiative transfer, even in vacant areas. This is an advantage of the present bidirectional upwind difference scheme. It also indicates that when we observe the outgoing upward radiance at a large zenith angle from above the cloud layer (such as μ = 0.55), the upward radiance leaving from the lateral side opposite to the solar incident side is comparable in intensity to that reflected on the cloud top. The qualitative pattern of the spatial distribution of the upward radiance within the cloud did not change with changes in the direction of the radiance. The maximum intensity appeared near the top of the cloud layer and at optically thick positions. On the other hand, the distribution pattern of the downward radiance differed with the different directions. The radiance at μ = −0.55, which is close to the solar incident angle, had a maximum in the skin layer of the cloud on the solar incident side and near the cloud top due to the strong forward scattering by the cloud particles.

Fig. 5.

The spatial distributions of the normalized intensities of the (top) upward diffusive radiance and (bottom) the downward diffusive radiance using the cloud model shown in Fig. 2a. The left-hand column is for μ = ±0.95 and the right-hand column is for μ = ±0.55. The cosine of the solar incident zenith angle is −0.6 and the azimuth angle is 0°.

Fig. 5.

The spatial distributions of the normalized intensities of the (top) upward diffusive radiance and (bottom) the downward diffusive radiance using the cloud model shown in Fig. 2a. The left-hand column is for μ = ±0.95 and the right-hand column is for μ = ±0.55. The cosine of the solar incident zenith angle is −0.6 and the azimuth angle is 0°.

Next, the calculation results for a cloud layer containing a more complicated spatial distribution are discussed. Figure 6 shows the distribution of the diffusive radiances at the upward directions (μ = 0.55) and downward directions (μ = −1.0) in the principal plane of ϕ = 0°, calculated for the cloud model shown in Fig. 2b. The geometrical thickness of the cloud layer was 600 m (the cloud-base height was 100 m, and the cloud-top height was 700 m). The maximum extinction coefficient of each vertical atmospheric column occurred at a height of 500 m, close to the top of the cloud layer. This feature is common to general stratocumulus cloud layers. It can be seen that the distribution of the upward radiance at μ = 0.55 was smoothed, and did not have a clear correspondence to the horizontal distribution of the cloud extinction coefficients. The distribution of the downward radiance at the nadir angle contained a larger variability than the upward radiance did, especially under the optically thin positions. However, the relationship between the optical thickness of the cloud and the downward radiances is not a simple one. This is a reason for the difficulties encountered in cloud remote sensing from observing the downward radiances.

Fig. 6.

The spatial distributions of the normalized intensities using the cloud model shown in Fig. 2b: (a) the upward diffusive radiance at μ = 0.55, and (b) the downward diffusive radiance at μ = −1.0. The cosine of the solar incident zenith angle is −0.6 and the azimuth angle is 0°.

Fig. 6.

The spatial distributions of the normalized intensities using the cloud model shown in Fig. 2b: (a) the upward diffusive radiance at μ = 0.55, and (b) the downward diffusive radiance at μ = −1.0. The cosine of the solar incident zenith angle is −0.6 and the azimuth angle is 0°.

We compared our method with data obtained using the Monte Carlo method. The spatial grid for both methods was set to 16.67 m along the x coordinate and 12.5 m along the z coordinate. The number of photons was determined to limit the computation errors to within 0.5% in the flux convergence. Furthermore, the phase function was not truncated in the Monte Carlo calculations. Thus, we consider the Monte Carlo results to be benchmark values in our comparison. Figure 7 shows the difference between the upward flux and the total downward flux calculated using our method and those obtained using the Monte Carlo method for the cloud model shown in Fig. 2a. The data in Fig. 7 show that the general features of the spatial distribution of the fluxes calculated using both methods were similar to each other. The upward flux and the total downward flux using our method are almost the same as the Monte Carlo results: the root-mean-square (RMS) differences divided by the corresponding mean value over the field are 0.024 for the upward flux and 0.017 for the total downward flux, respectively. Figure 8 illustrates the difference between the intensities of upward radiance (μ = 0.95) and the downward radiance (μ = −0.95) with the azimuth angle of ϕ = 0°. This indicates that the intensities of radiance by using our method are also in good agreement to those calculated by the Monte Carlo method. However, the discrepancy in the downward radiance is relatively large, especially near the sunny side of the cloud laterals, where the angular distribution of the radiation was highly anisotropic because of less scattering by the highly asymmetric phase function of cloud particles. The root-mean-square differences divided by the corresponding mean value over the field are 0.063 for the upward radiance and 0.079 for the downward radiance, respectively. These comparisons to the Monte Carlo method prove that our method is sufficient to simulate radiation fields of both fluxes and radiances in inhomogeneous atmospheres.

Fig. 7.

The spatial distributions of the differences using the Monte Carlo data from (a) the normalized upward flux and (b) the normalized total downward flux for the cloud model shown in Fig. 2a. The negative values (underestimation using the present method) are shown by the dotted line. The spatial resolution for both methods is dx = 16.67 m and dz = 12.5 m.

Fig. 7.

The spatial distributions of the differences using the Monte Carlo data from (a) the normalized upward flux and (b) the normalized total downward flux for the cloud model shown in Fig. 2a. The negative values (underestimation using the present method) are shown by the dotted line. The spatial resolution for both methods is dx = 16.67 m and dz = 12.5 m.

Fig. 8.

The spatial distributions of the differences using the Monte Carlo data from (a) the normalized upward radiance (μ = 0.95) and (b) the normalized downward radiance (μ = −0.95) for the cloud model shown in Fig. 2a. The negative values (underestimation using the present method) are shown by the dotted line. The cosine of the solar incident zenith angle is −0.6 (θ0 = −53°) and the azimuth angle of radiances is 0°. The spatial resolution for both methods is dx = 16.67 m and dz = 12.5 m.

Fig. 8.

The spatial distributions of the differences using the Monte Carlo data from (a) the normalized upward radiance (μ = 0.95) and (b) the normalized downward radiance (μ = −0.95) for the cloud model shown in Fig. 2a. The negative values (underestimation using the present method) are shown by the dotted line. The cosine of the solar incident zenith angle is −0.6 (θ0 = −53°) and the azimuth angle of radiances is 0°. The spatial resolution for both methods is dx = 16.67 m and dz = 12.5 m.

## 4. Resolution and accuracy

### a. Spatial resolution

In general, the time required for computing a solution increases with the increasing accuracy required, and it depends on the spatial and angular resolutions of the discretization. For efficient calculations, it is necessary to know how the calculation results are influenced by the resolution. We first have to calculate the radiation fields for the cloud field shown in Fig. 2a by changing the spatial resolution. Table 2 shows the RMS difference divided by the mean value over the entire domain in the upward flux, total downward flux, net flux along the x coordinate, and the radiances at the μ = 0.95 (upward) and μ = −0.95 (downward) directions with an azimuthal angle of ϕ = 0°. For the reference values, we selected the solution computed with an x-coordinate interval (dx) of 16.67 m and a z-coordinate interval (dz) of 12.5 m. The maximum degree (L) and order (M) of the spherical harmonic series expansion were set to values of 15 and 8, respectively. The normalized RMS difference between the Monte Carlo results and the reference values is also listed in Table 2. The results suggest that in this example, when the x-coordinate resolution is halved, the normalized RMS differences of the upward, the downward, and the net flux values are also nearly halved, except for the upward flux with dz = 12.5 m. On the other hand, the z-coordinate resolution is less significant for improvement in accuracy of the flux values than the x-coordinate resolution. Accuracy of the intensities of upward and downward radiance is dependent on both the x- and z-coordinate resolution. The large discrepancies, especially in the radiance comparisons, occur locally on the sunny side of the cloud laterals and the downstream positions of the direct solar incident radiance that grazes the cloud edges. In addition, the number of iterations (389 for the reference resolution case) increases with increasing number of grid points. This means that when the spatial resolution becomes finer, more calculation time is required because of the increase in the number of grid points and an increase in the number of iterations needed for convergence. Figure 9 shows the reduction in the |N(r)n+1N(r)n|/|N(r)n| values with increasing number of iterations using the SOR method in the case where dx = 50 m and dz = 25 m, and in the case where dx = 25 m and dz = 12.5 m. This suggests that at the beginning of the iterations, the reduction rate is almost independent of the spatial resolution, but the rate in the case where dx = 25 m and dz = 12.5 m decreases with increasing number of iterations compared to the case where dx = 50 m and dz = 25 m.

Table 2.

The RMS difference divided by the mean value over the domain in the upward flux, the total downward flux, the net flux along the x coordinate, and the upward and downward radiances with a given spatial resolution. The reference values are for dx = 16.67 m and dz = 12.5 m. The maximum degree (L) and order (M) of the spherical harmonic series expansion were set to values of 15 and 8, respectively. The azimuth angle of the radiances was set to 0°. The RMS difference between the MC result and the reference is also listed.

Fig. 9.

Values of |N(r)n+1N(r)n|/|N(r)n| as a function of the number of iterations using the SOR method. The cloud model used is shown in Fig. 2a. The solid line shows the case where dx = 50 m and dz = 25 m. The dotted line shows the case where dx = 25 m and dz = 12.5 m in Table 2.

Fig. 9.

Values of |N(r)n+1N(r)n|/|N(r)n| as a function of the number of iterations using the SOR method. The cloud model used is shown in Fig. 2a. The solid line shows the case where dx = 50 m and dz = 25 m. The dotted line shows the case where dx = 25 m and dz = 12.5 m in Table 2.

### b. Angular resolution

In this section, we will discuss the dependence of the accuracy on the angular resolution. Figure 10 shows the RMS differences divided by the corresponding mean values over the entire domain in the upward flux and total downward flux of the cloud model shown in Fig. 2a, as a function of L (the maximum number of degrees) and M (the maximum number of orders) of the spherical harmonic series. Figure 11 is the same but for the radiance at the μ = 0.95 (upward) and μ = −0.95 (downward) directions with an azimuthal angle of ϕ = 0°. The solution computed with L = 23 and M = 14 was selected for the reference values. The reference value was closest to the Monte Carlo data in the results computed with smaller L or M values. The data in Fig. 10 suggest that even using a smallest number, such as L = 11 and M = 4, is sufficient for the convergence of the upward and total downward flux to within 1%. The RMS difference in the upward flux effectively decreases with increasing M, and dependence on L is small. In contrast, the accuracy of the total downward flux is dependent on both L and M. On the other hand, Fig. 11 shows that the normalized RMS difference in the radiation intensities increases more significantly with decreasing values of L, whereas M is less significant for accuracy of calculation. In addition, the number of iterations for convergence of the solution was almost independent of the values of L and M. The solution for the cloud model shown in Fig. 2a was obtained after about 180 iterations for all cases of L and M. This means that the increase in the calculation time for a finer angular resolution was only due to the increase in the vector and matrix size.

Fig. 10.

Difference in the normalized root mean squares in (left) the upward flux and (right) the total downward flux from those calculated using L = 23 and M = 14, as a function of L and M for the spherical harmonic series. The cloud model used for calculation is shown in Fig. 2a.

Fig. 10.

Difference in the normalized root mean squares in (left) the upward flux and (right) the total downward flux from those calculated using L = 23 and M = 14, as a function of L and M for the spherical harmonic series. The cloud model used for calculation is shown in Fig. 2a.

Fig. 11.

Difference in the normalized RMSs in the intensity of radiance for (left) μ = 0.95 and (right) μ = −0.95 from those calculated using L = 23 and M = 14, as a function of L and M for the spherical harmonic series. The cloud model used for calculation is shown in Fig. 2a and the azimuth angle of radiances is 0°.

Fig. 11.

Difference in the normalized RMSs in the intensity of radiance for (left) μ = 0.95 and (right) μ = −0.95 from those calculated using L = 23 and M = 14, as a function of L and M for the spherical harmonic series. The cloud model used for calculation is shown in Fig. 2a and the azimuth angle of radiances is 0°.

## 5. Summary and conclusions

In this study, we have developed a new calculation scheme for the explicitly discretized solution of the three-dimensional radiation transfer equation in an inhomogeneous atmosphere. Discretization techniques are used to separate the independent variables involved in the 3D RTE. The solution employed a spherical harmonic series expansion to discretize the terms that depend on the direction of the radiance. The spherical harmonic transform is appropriate for obtaining the intensity of the radiation in an arbitrary direction. We have also applied the finite-volume method to discretize the terms that depend on the spatial coordinates. For the discretization of the partial-differential terms, we developed a bidirectional upwind difference scheme, which is a specialized scheme for the RTE expanded by the spherical harmonic functions. This scheme treats radiation fields by separating each direction. This scheme can make the discretized 3D RTE determinate and a convergent solution can be obtained. The 3D RTE can be formulated as a simultaneous linear equation, which can be described in a vector–matrix form with a sparse matrix. We use the successive overrelaxation (SOR) method to solve the discretized equation. Our solution is useful and effective for stable and efficient computation of the radiation fields in inhomogeneous clouds.

To show the good performance of our scheme, we have shown several examples of the radiation field calculated for two-dimensional inhomogeneous cloud models. We have shown that our method can properly and efficiently simulate the radiation field of both the flux and radiance in inhomogeneous clouds. Validation of our calculation results was carried out by comparing our data with data obtained using the Monte Carlo method. The upward flux and the total downward flux (direct plus diffusive) obtained using our method were in good agreement with data obtained using the Monte Carlo method, to within a few percent of the normalized root-mean-square difference. However, the downward radiance exhibited rather large discrepancies at the positions where the angular distribution of the radiation was anisotropic due to less scattering by the highly asymmetric phase function.

We calculated the flux and radiance for various spatial and angular resolutions to estimate the errors caused by our discretizations. Our calculations indicate that the horizontal resolution is more responsible for the errors in the fluxes than the vertical resolution, whereas both resolutions are effective to the errors in the upward and downward radiances. Accuracy in the upward and the total downward fluxes is weakly dependent on the angular resolution. However, the radiances are quite influential on the angular resolution, especially at the maximum degree (L) of the spherical harmonic series.

We believe that the present solution is efficient and computationally stable. Our method will be helpful for investigating and quantitatively estimating the effect of cloud spatial inhomogeneity on the corresponding radiation field. Such studies are required to improve cloud simulations in high-resolution models, as well as to develop remote sensing techniques for horizontally inhomogeneous clouds and the earth’s surface. Our method can be further refined. Other types of spatial discretization, such as hexagonal or tetrahedral, can be applied to the finite-volume method. These techniques may be more appropriate in representing a cloudy atmosphere and the corresponding radiation fields. The upwind difference scheme used in this study is a first-order precision scheme. A more refined upwind difference scheme may enable us to obtain a more accurate solution. We intend to carry out improvements to our present solution in future work.

## Acknowledgments

We thank Dr. H. Iwabuchi for providing the data of the Monte Carlo simulations. This study was partly supported by the Ministry of Education, Culture, Sports, Science and Technology, under a Grant-in-Aid for Scientific Research [(A)17204039].

## REFERENCES

REFERENCES
Barker
,
H. W.
, and
J. A.
Davies
,
1992
:
Cumulus cloud radiative properties and the characteristics of satellite radiance wavenumber spectra.
Remote Sens. Environ.
,
42
,
51
64
.
Cahalan
,
R. F.
, and
J. B.
Snider
,
1989
:
Marine stratocumulus structure.
Remote Sens. Environ.
,
28
,
95
107
.
Cahalan
,
R. F.
,
W.
Ridgway
,
W. J.
Wiscombe
, and
S.
Gollmer
, and
Harshvardhan
,
1994
:
Independent pixel and Monte Carlo estimates of stratocumulus albedo.
J. Atmos. Sci.
,
51
,
3776
3790
.
Cahalan
,
R. F.
, and
Coauthors
,
2005
:
The 13RC: Bringing together the most advanced radiative transfer tools for cloudy atmospheres.
Bull. Amer. Meteor. Soc.
,
86
,
1275
1293
.
Chai
,
J. C.
,
P-f
Hsu
, and
Y. C.
Lam
,
2004
:
Three-dimensional transient radiative transfer modeling using the finite-volume method.
,
86
,
299
313
.
Evans
,
K. F.
,
1993
:
Two-dimensional radiative transfer in cloudy atmospheres: The spherical harmonic spatial grid method.
J. Atmos. Sci.
,
50
,
3111
3124
.
Evans
,
K. F.
,
1998
:
The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer.
J. Atmos. Sci.
,
55
,
429
446
.
Gabriel
,
P. M.
,
S.
Tsay
, and
G. L.
Stephens
,
1993
:
A Fourier–Riccati approach to radiative transfer. Part I: Foundations.
J. Atmos. Sci.
,
48
,
2436
2447
.
Gu
,
Y.
, and
K. N.
Liou
,
2001
:
Radiation parameterization for three-dimensional inhomogeneous cirrus clouds: Application to climate models.
J. Climate
,
14
,
2443
2457
.
Hartman
,
C. M.
, and
J. Y.
Harrington
,
2005
:
Radiative impacts on the growth of drops within simulated marine stratocumulus. Part I: Maximum solar heating.
J. Atmos. Sci.
,
62
,
2323
2338
.
Ishida
,
H.
, and
S.
Asano
,
2007
:
A quasi-analytic solution of the radiative transfer equation for three-dimensional atmospheres.
,
103
,
371
393
.
Iwabuchi
,
H.
, and
T.
Hayasaka
,
2002
:
Effects of cloud horizontal inhomogeneity on the optical thickness retrieved from moderate-resolution satellite data.
J. Atmos. Sci.
,
59
,
2227
2242
.
Kobayashi
,
T.
,
1991
:
Reflected solar flux for horizontally inhomogeneous atmospheres.
J. Atmos. Sci.
,
48
,
2436
2447
.
Marshak
,
A.
, and
A. B.
Davis
,
2005
:
3D Radiative Transfer in Cloudy Atmospheres.
Springer, 686 pp
.
Moeng
,
C.
,
P. P.
Sullivan
, and
B.
Stevens
,
1999
:
Including radiative effects in an entrainment rate formula for buoyancy-driven PBLs.
J. Atmos. Sci.
,
56
,
1031
1049
.
Ramanathan
,
V.
,
R. D.
Cess
,
E. F.
Harrison
,
P.
Minnis
,
B. R.
Barkstrom
,
E.
, and
D.
Hartmann
,
1989
:
Science
,
243
,
57
63
.
Stephens
,
G. L.
,
1988a
:
Radiative transfer through arbitrarily shaped optical media. Part I: A general method of solution.
J. Atmos. Sci.
,
45
,
1818
1836
.
Stephens
,
G. L.
,
1988b
:
Radiative transfer through arbitrarily shaped optical media. Part II: Group theory and simple closures.
J. Atmos. Sci.
,
45
,
1837
1848
.
Varnai
,
T.
, and
A.
Marshak
,
2002
:
Observations of three-dimensional radiative effects that influence MODIS cloud optical thickness retrievals.
J. Atmos. Sci.
,
59
,
1607
1618
.
Wiscombe
,
W. J.
,
1977
:
The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions.
J. Atmos. Sci.
,
34
,
1408
1422
.
Wu
,
X.
, and
X-Z.
Liang
,
2005
:
Radiative effects of cloud horizontal inhomogeneity and vertical overlap identified from a monthlong cloud-resolving model simulation.
J. Atmos. Sci.
,
62
,
4105
4112
.

### APPENDIX A

#### Definition of the Coefficients for the Hemispherical Integral

The components of the spherical harmonic-transformed radiance for a given direction, such as J+e;c/s,l,m or Je;c/s,l,m in Eq. (13), are defined by a hemispherical integral of the radiances using the spherical harmonic function. For example, the transformed radiance along the positive x-coordinate direction (i.e., the right-hand direction in Fig. 1b), J+e;c/s,l,m, is defined by

which means the integral is in the domain of [−1, 1] for μ and [−π/2, π/2] for ϕ. Substituting N(μ, ϕ) from Eq. (2) into Eq. (A1), means that Eq. (A1) can be expressed by the coefficients for the hemispherical integral, Rc/s;l,m;l′,mrgt in the form

where

to obtain Eq. (14). In a similar way, the component of the transformed radiance for the negative direction along the x coordinate (i.e., left-hand direction of Fig. 1b), Je;c/s,l,m, is defined by the hemispherical integral in the domain of [−1, 1] for μ and [π/2, 3π/2] for ϕ as

where

to obtain Eq. (15). From Eqs. (A2) and (A4), J+e;c/s,l,m and Je;c/s,l,m satisfy the relationship

The coefficients of the hemispherical integral for the positive direction along the y coordinate are also defined by the hemispherical integral in the domain of [−1, 1] for μ and [0, π] for ϕ as

and those for the negative direction are defined by the integral in the domain of [−1, 1] for μ and [π, 2π] for ϕ as

The coefficients of the hemispherical integral for the positive and negative directions along the z coordinate are defined by the hemispherical integrals in the domain of [0, 1] for μ and [0, 2π] for ϕ, respectively, as

and in the domain [−1, 0] for μ and [0, 2π] for ϕ as

### APPENDIX B

#### Reflecting Surface Boundary Conditions

When radiance at the bottom boundary is N(μ, ϕ) and the bidirectional reflection function of the surface is given by Rs(μ, ϕ; μ′, ϕ′), the components of the spherical harmonic-transformed radiance reflected by the surface, Nr;c/s,l,m, is written as

where

Equation (B1) can be written in the vector–matrix form

where vector Nr consists of Nr;c/s,l,m and the elements of matrix 𝗥s are Rc/s;l,m;l′,ms.

## Footnotes

Corresponding author address: H. Ishida, Research and Information Center, Tokai University, 2-28-4, Tomigaya, Shibuya-ku, Tokyo 151-0063, Japan. Email: ishidah@yoyogi.ycc.u-tokai.ac.jp