## 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

Many satellite observations [e.g., the pioneering study of Ramanathan et al. (1989)] have shown that clouds have a large effect on the radiation budget of the atmosphere. Numerical simulations also suggest that radiation processes are very important for cloud formation and decay (e.g., Moeng et al. 1999; Hartman and Harrington 2005). However, radiative transfer within clouds is very complicated and difficult to evaluate because of spatial inhomogeneities, such as the existence of complicated geometric shapes, the large variability in the extinction coefficient within a cloud, and the variety of water phases present. Inhomogeneity can also be present in a natural cloud layer that seems to be horizontally uniform in its appearance, such as marine stratocumulus (e.g., Cahalan and Snider 1989). Past simulation studies have shown that involving the effect of a spatially inhomogeneous cloud leads to radiative transfer calculation data that are different from those using a simple plane-parallel approximation (e.g., Barker and Davies 1992; Cahalan et al. 1994). Omitting the horizontal inhomogeneity of cloud layers in a general circulation model (GCM) results in an incorrect shortwave radiation flux (e.g., Wu and Liang 2005). In addition, recent advances in atmospheric remote sensing devices (either satellite or ground based) enable radiance measurements to be performed with a finer spatial resolution than previously, but observations using such devices are greatly influenced by the spatial structures of clouds (e.g., Varnai and Marshak 2002). Therefore, an effective method for calculating radiative transfer in spatially inhomogeneous atmospheres is required to investigate the radiative properties of actual clouds, to provide parameterizations of the radiative effect of inhomogeneous clouds in cloud formation models and GCMs, and to improve the accuracy of remote sensing products.

In general, the radiative transfer equation (RTE) in a three-dimensional (3D) medium is a five-dimensional integral-differential equation (i.e., three variables for the spatial coordinates and two angular variables for the direction of radiance), which results in a boundary value problem. The incident radiances on the boundaries are given, and the outgoing radiances need to be calculated. Various methods for calculating the radiation field in a 3D atmosphere have been proposed because it is difficult to solve a 3D RTE analytically, unlike the plane-parallel (one-dimensional) approximation. Cahalan et al. (2005) reviewed recent results and compared the different performances. The Monte Carlo (MC) method, which statistically simulates the transport behavior of photons in an atmosphere by tracking photon paths, is a highly effective and appropriate method for calculating 3D radiation fields (e.g., Iwabuchi and Hayasaka 2002). Another type of solution using analytic (or deterministic) methods has also been developed. In contrast to the MC method, analytic methods are able to explicitly represent the RTE and can calculate the entire radiation field for a given domain. An advantage of analytic methods is that they can explicitly represent the relationship between the horizontal structure of a cloud layer and the corresponding radiation field. Therefore, analytic methods can be used to investigate the effect of cloud inhomogeneity on the radiation transfer process. For example, applying the solution method of Stephens (1988a), Stephens (1988b) introduced a formulation to distinguish various groups of spatial scales of variations that make up the radiance field.

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 *N*^{c/s}_{l,m} is the transformed intensity with spherical harmonics functions, *Y*^{c/s}_{l,m}(*μ*, *ϕ*), that are defined by

where *P*^{m}_{l}(*μ*) 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*.

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 *Y*^{c/s}_{l,m}(*μ*, *ϕ*) results in

where

and the *a _{l}*

_{′}

_{,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

*a*

_{l}_{′}

_{,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, *P _{l}*(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 *a*_{l,m}, *b*_{l,m}^{c/s}, and *c*_{l,m}, respectively, and the elements of matrix 𝗖 are *α* + *s**χ*_{l}/(2*l* + 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, *N*^{c/s}_{l,m}, on each boundary of the cell (e, w, b, t in Fig. 1b). The *N*^{c/s}_{l,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 *z*–*y* plane (corresponding to the plane including position e in Fig. 1b), *N*^{c/s}_{l,m} can be written as

where *J*^{+}_{e;c/s,l,m} and *J*^{−}_{e;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 *J*^{−}_{e;c/s,l,m}, can be expressed by

and

where *R*_{c/s;l,m;l′,m′}^{rgt} and *R*_{c/s;l,m;l′,m′}^{lft} 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.

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 **r**_{P} is the position vector at grid point P. On the other hand, *J*^{−}_{e;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 **r**_{E} 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 **J**^{−}_{w}, are given by

and

where **r**_{W} 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}, **J**^{−}_{n}, **J**^{+}_{s}, and **J**^{−}_{s}, are given by

where **r**_{N} and **r**_{S} 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}, **J**^{−}_{t}, **J**^{+}_{b}, and **J**^{−}_{b}, are given by

where **r**_{T} and **r**_{B} are the position vectors at grid points T and B, respectively. The elements of the matrices, *R*_{c/s;l,m;l′,m′}^{frn}, *R*_{c/s;l,m;l′,m′}^{bck}, *R*_{c/s;l,m;l′,m}^{up}, and *R*_{c/s;l,m;l′,m}^{dwn}, are also given in appendix A.

Next, we take into account the conservation of a physical quantity (corresponding here to *N*^{c/s}_{l,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 P_{b} 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 P_{c} 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+1} − *N*(*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 2*n*^{2} + {2(*L* + 1)(*M* + 2) + 2(*L* + 3) + *L* + 3 + (*L* + 1)(*M*/2 + 1) + 3} × *n*, where *n* = (*L* + 1) × (*M* + 1) × *n _{x}* ×

*n*×

_{y}*n*and

_{z}*n*×

_{x}*n*×

_{y}*n*is the number of spatial grid points. After

_{z}**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 *x*–*z* 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.

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.

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.

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.

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.

## 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 d*z* = 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+1} − *N*(*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.

### 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.

## 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

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### 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 *J*^{−}_{e;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, *R*_{c/s;l,m;l′,m′}^{rgt} 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), *J*^{−}_{e;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 *J*^{−}_{e;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 *R ^{s}*(

*μ*,

*ϕ*;

*μ*′,

*ϕ*′), the components of the spherical harmonic-transformed radiance reflected by the surface,

*N*

_{r}_{;}

_{c}_{/}

*,*

_{s}*,*

_{l}*, is written as*

_{m}where

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

where vector **N*** _{r}* consists of

*N*

_{r}_{;}

_{c}_{/}

*,*

_{s}*,*

_{l}*and the elements of matrix 𝗥*

_{m}^{s}are

*R*

_{c/s;l,m;l′,m′}

^{s}.

## 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