Abstract

Cloud-resolving models—in particular, large-eddy simulation (LES) models—are important tools to improve the understanding of cloud–radiation interactions. A method is presented for accurate, yet fast, three-dimensional calculation of surface shortwave irradiance within an LES model using the tilted independent column approximation with smoothing of the diffuse irradiance. The algorithm calculates a tilted optical thickness for each surface pixel that is then used as input to a one-dimensional radiative transfer code. In a sensitivity analysis, it is shown that this calculation can even be replaced by a simple precalculated lookup table that tabulates surface irradiance as a function of only solar zenith angle and cloud optical thickness. Because the vertical variability of the cloud is of little relevance for the surface irradiance, this approximation introduces little extra uncertainty. In a final step, surface irradiance is smoothed to account for horizontal photon transport between individual columns. The algorithm has been optimized for parallelization, which enhances its applicability in LES models. In this implementation, the total computational time of the LES model increased by only 3% relative to the reference run without radiation. Comparisons between the fast approximation and detailed three-dimensional radiative transfer calculations showed very good agreement for different cloud conditions and several solar zenith and azimuth angles, with a root-mean-square difference of 6%.

1. Introduction

Clouds generally increase the atmospheric albedo and reduce the amount of solar radiation reaching the surface. This leads to a reduction of the upward heat flux from the surface into the atmosphere, which has in particular an effect on boundary layer clouds: changes in buoyancy within the boundary layer, which drive the turbulent convective motions, feed back on the formation of clouds. As shown by Schumann et al. (2002), the consideration of three-dimensional (3D) radiative effects (e.g., surface heating rates) can have significant influence on the development of convective clouds. Radiation may affect clouds through different mechanisms. In the longwave spectral range, two effects happen: emission of longwave radiation from cloud top cools while the cloud-base absorption of longwave emission from the surface leads to a warming. We concentrate on the influence of differential surface heating on clouds caused by the inhomogeneous distribution of irradiance at the ground—in particular that due to cloud shadows.

The influence of differential surface heating due to anvil shadow underneath a supercell on its structure and dynamics has been investigated by Markowski and Harrington (2005). Markowski and Harrington (2005) as well as Schumann et al. (2002) used coarse approximations to account for differential surface heating. To investigate these effects further, an accurate yet fast method to implement differential surface irradiation in a cloud-resolving model (CRM) is presented.

The radiation reaching the surface is composed of transmitted shortwave radiation and longwave radiation emitted by the atmosphere. The relevance of these two contributions is investigated. For this purpose, shortwave and longwave surface irradiance underneath a broken cloud field, calculated with the large-eddy simulation model known as “EULAG” (from Eulerian–Lagrangian), has been calculated with the three-dimensional radiative transfer code known as the Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres (MYSTIC; Mayer 2000a, b). The broken cloud field has a cloud cover of 0.3, the cloud base is at 600 m, and the cloud top is at 1600 m. The cloud field and the results of the radiation calculations are shown in Fig. 1. The following can be seen: the distribution of longwave irradiance is much smoother than that of the shortwave irradiance. This situation is due to the longwave radiation being emitted in all directions whereas shortwave irradiance is dominated by the direct solar radiation. The difference between the lowest and highest value of the solar irradiance is approximately 1000 W m−2; that of the thermal irradiance is approximately 30 W m2. Thus the variability of longwave irradiance amounts to only 3% of that in the solar spectral range. For the purpose of studying the impact of inhomogeneous surface irradiance on cloud formation, we can therefore safely restrict ourselves to the shortwave spectral range.

Fig. 1.

Integrated (right) longwave and (middle) shortwave irradiance underneath a (left) cloud field (the labeling of the axes corresponds to the model grid) with a solar zenith angle of 45° (the sun is in the south; at the bottom of the figure); note the different scales of the two rightmost images.

Fig. 1.

Integrated (right) longwave and (middle) shortwave irradiance underneath a (left) cloud field (the labeling of the axes corresponds to the model grid) with a solar zenith angle of 45° (the sun is in the south; at the bottom of the figure); note the different scales of the two rightmost images.

Surface irradiance depends on the geometry of the clouds, the spatial distribution of the liquid water content, and the droplet size distribution within the clouds (e.g., O’Hirok and Gautier 1998). Typical three-dimensional phenomena are the following: 1) in the shadow of the clouds, the amount of solar radiation reaching the surface is considerably reduced relative to in nonshadow regions—in particular, for clouds with large water content; 2) in nonshadow regions, more solar radiation is reaching the surface relative to in the cloudless case because of scattering from illuminated cloud sides; and 3) the shadow occurs not directly below the cloud but is shifted to the side depending on the solar zenith and azimuth angle. An accurate algorithm should be able to represent these phenomena. Three-dimensional radiative transfer codes represent an accurate way to calculate the irradiance distributions, but they have high computational demands (Cahalan et al. 2005), a fact that usually prevents their direct application in CRMs. For this purpose, fast parameterizations are required that nevertheless correctly include three-dimensional effects.

A widely employed method is the independent column approximation (ICA). In particular in the field of remote sensing, it is often denoted by the term “independent pixel approximation” (Cahlan et al. 1994). The atmosphere is divided into single columns, and for each the radiative transfer calculations are performed independently. The ICA can yield reasonable results for mean values over large domains but not for the calculation of horizontal distributions. The neglect of horizontal transport between columns causes errors. In contrast to 3D radiative transfer calculations, the cloud shadows calculated with ICA are always located directly underneath the cloud irrespective of the solar zenith angle, which is obviously wrong for any zenith angle except zero. In addition to the shifted location, the ICA shadow is systematically too small. In the real world, the area of the shadows increases with increasing solar zenith angle under a broken cloud field, whereas in the ICA the shadow area does not depend on solar zenith angle. If one considers as a simple and extreme example a chessboard field of cloud cubes (see, e.g., Weinman and Harshvardhan 1982), no direct radiation would reach the ground for solar zenith angles larger than 45°, whereas in the ICA one-half of the area would be sunlit for any solar zenith angle. This is illustrated in the top panels of Fig. 2. The overestimation of the direct irradiance in clear sky is connected to an underestimation of the diffuse irradiance: because less radiation is scattered out of the direct beam in the ICA, the corresponding diffuse irradiance must be smaller. The two effects partially cancel, but the spatial distribution of the total irradiance at the ground is obviously completely different from reality. The spatial distribution of the irradiance underneath the chessboard cloud field is shown in the bottom panels of Fig. 2 for a sun zenith angle of 45°. In the ICA calculations, direct radiation reaches the surface by passing between the clouds and the shadows appear at the same position as the clouds (e.g., in the lower-left corner). In the three-dimensional calculations, direct radiation reaches the surface only in small areas where it penetrates through the corners of the clouds (see dashed lines in the top-left panel of Fig. 2).

Fig. 2.

(top) Schematic of cloud shadows for an idealized cloud field with a chessboard pattern for a solar zenith angle of 45°: (left) 3D and (right) ICA. (bottom) Total surface irradiance underneath the chessboard cloud field: (left) 3D and (right) ICA calculations for a solar zenith angle of 45° (the sun is in the south).

Fig. 2.

(top) Schematic of cloud shadows for an idealized cloud field with a chessboard pattern for a solar zenith angle of 45°: (left) 3D and (right) ICA. (bottom) Total surface irradiance underneath the chessboard cloud field: (left) 3D and (right) ICA calculations for a solar zenith angle of 45° (the sun is in the south).

The ICA is therefore only suitable for homogeneous overcast clouds or for low-resolution models for which the assumption of independent columns is applicable, that is, for which the grid boxes are large enough so that the shadow falls in the same grid box anyway. In CRMs, the ICA certainly cannot correctly represent cloud shadows for nonzero solar zenith angles. In this paper, a new method for the accurate yet fast calculation of surface irradiance within a CRM is presented that is based on the tilted independent column approximation (TICA). The “tilted independent pixel approximation” is a one-dimensional (1D) radiative transfer approximation introduced by Várnai and Davies (1999), who examined the processes through which cloud heterogeneities influence solar radiation and explained the relevant mechanisms. Here we use the term tilted independent column (rather than pixel) approximation, which is more appropriate for application in three-dimensional grids. The TICA treats the direct beam exactly while approximating the diffuse radiation. The method guarantees that at least the shadows are at the correct location and, as we will show, may also give accurate values for the diffuse component if the diffuse irradiance is smoothed over an appropriate distance. Here we describe a very efficient implementation of the TICA that can be implemented in CRMs at virtually no cost.

The paper is organized as follows. Section 2 explains the TICA and the realization within parallelized CRMs. In section 3, the TICA is validated by comparison with detailed three-dimensional radiative transfer calculations. Section 4 gives a summary and some concluding remarks. The effects on cloud formation and evolution will be presented in a companion paper.

2. Method

Our main motivation for developing the fast radiation algorithm was to provide an accurate radiation parameterization for the anelastic nonhydrostatic large-eddy simulation model EULAG (Smolarkiewicz and Margolin 1997; Grabowski and Smolarkiewicz 2002). In its standard version, no radiation is included. The heat flux at the surface is usually specified as a constant value over the model domain. We included the new algorithm, based on the TICA, in EULAG to calculate surface solar irradiance as accurately, yet as quickly, as possible. The feedback of surface solar irradiance on cloud formation is then simulated by adaptation of the surface heat flux to the incident solar irradiance, assuming that the surface cools quickly when the sun is masked by clouds. The 3D cloud fields used for Fig. 1 and for the validation of the new algorithm (shown later) are generated with EULAG. These simulations were performed with a horizontal resolution of 50 m and a vertical resolution of 40 m within a model domain of 6.4 km × 6.4 km. In the following sections, we describe the TICA and its practical implementation and introduce some improvements in terms of speed and accuracy.

a. Tilted independent column approximation

In the TICA, the radiative transfer calculations are performed independently for single columns that are tilted in the direction of the sun, along the solar zenith and azimuth angles. A schematic of the TICA method is shown in Fig. 3 in which a thick line underneath the cloud shows the position of the cloud shadow. The left panel of Fig. 3 illustrates some real photon paths. Photons can reach the surface directly, be reflected by cloud sides, or be scattered within clouds and thereafter reach the surface or escape at the top of the atmosphere. The middle panel of Fig. 3 illustrates how these photons would be treated in a conceptual version of the TICA: periodic boundary conditions are assumed within each column, which prevent photons from leaving each column. (Note that this is only for illustration; in the real implementation, each column is used independently as input to a 1D radiative transfer code.) The direct photon path is the same as in the three-dimensional treatment, whereas the scattered photons are confined to the column. However, the location and size of the shadow is the same in both cases. Although direct irradiance is calculated correctly in the TICA, the calculation within independent columns leads to incorrect diffuse irradiance. In reality, different effects occur; for example, the diffuse irradiance outside the shadows is often larger than the corresponding cloudless sky value (e.g., Mims and Frederick 1994) because of reflection of radiation by cloud sides. In contrast, the diffuse irradiance assumes its cloudless sky value in the TICA outside cloud shadows. The missing horizontal photon transport between columns in the TICA causes too-high values of diffuse irradiance in cloudy columns and too-low values in cloudless columns. This error can be reduced considerably (as shown later) by applying a smoothing function to the calculated distribution of diffuse irradiance, yielding an approximated allowance of horizontal photon transport.

Fig. 3.

Schematic of (left) 3D calculation, (middle) TICA, and (right) ICA: possible paths are plotted.

Fig. 3.

Schematic of (left) 3D calculation, (middle) TICA, and (right) ICA: possible paths are plotted.

The actual technical realization of the TICA is illustrated in Fig. 4. To convert from the 3D cloud information to columns that can be used as input to a standard 1D radiative transfer code, the distribution of optical properties along individual rays is translated to a vertical column where the 3D grid boxes are translated to layers, the thickness of which reflects the pathlength in each 3D grid box. Using this atmosphere as input to a 1D model, the direct irradiance is calculated exactly but the diffuse irradiance is approximated, because horizontally infinite layers are assumed.

Fig. 4.

Calculation of the vertical profile of optical properties (gray tone) within a tilted column: (left) geometrical photon path and (right) the corresponding 1D radiative transfer model input.

Fig. 4.

Calculation of the vertical profile of optical properties (gray tone) within a tilted column: (left) geometrical photon path and (right) the corresponding 1D radiative transfer model input.

b. Smoothing of the diffuse component

The TICA treats the direct beam exactly while the diffuse radiation is approximated, because photons cannot move between adjacent tilted columns. This constraint then causes errors in the diffuse irradiance at the surface because each surface element receives diffuse radiation from all directions. A straightforward way to consider this is to average the diffuse irradiance over some smoothing kernel that depends on the clouds—in particular, the cloud-bottom height. If one assumes isotropic radiance, a surface element receives 50% of its irradiance from a cone of half-opening angle 45°. With the assumption of a cloud base at 1 km, this translates to a horizontal distance at cloud bottom of 2 km. A typical smoothing kernel would therefore have a width of about 2 km. This is actually larger than the individual clouds in our simulations, and we found that smoothing the diffuse irradiance over the complete model domain was actually as good as any other tested smoothing kernel. For the application in the CRM we therefore decided simply to average the diffuse irradiance over the complete model domain. For larger domain sizes of more irregular cloud fields, a more general approach using a smoothing kernel, such as the nonlocal independent pixel approximation by Marshak et al. (1998), would be required, which, however, would also be straightforward to implement into our model framework. The validity of our assumption will later be shown by comparison between our implementation of the TICA and a real 3D model.

c. Relevance of the cloud vertical distribution

To speed up the calculations even further, we studied the impact of the vertical distribution of cloud optical properties on surface irradiance with a 1D model, with the aim to replace the 1D radiative transfer calculations by a simple lookup table providing surface irradiance for given cloud optical thickness and solar zenith angle. For this study, profiles of pressure, temperature, water vapor, and so on were taken from the midlatitude summer profile by Anderson et al. (1986). Four different cloud scenarios have been examined (schematics are shown in Fig. 5). The first scenario studied the influence of the vertical location of a given cloud. To accomplish this, the height of a sample cloud has been varied while keeping its optical properties constant: One-dimensional calculations have been performed for a cloud with a geometrical thickness of 200 m and an optical thickness of 2 at a height of 200 m, 500 m, 1 km, 2 km, 5 km, and 7 km, respectively. This is illustrated in Fig. 5a, showing the same sample 1D cloud profile at different heights. The second scenario (Fig. 5b) examined the influence of the geometrical thickness of the cloud for a given optical thickness. For that purpose, the geometrical thickness has been varied for a cloud with constant cloud-base height at 1 km and a constant optical thickness of 2 (Fig. 5b). The cloud had a geometrical thickness of 200 m, 500 m, 1 km, 2 km, 3 km, and 5 km for the different calculations. The third scenario (Fig. 5c) examined the influence of the vertical distribution of liquid water within a cloud with constant base, geometrical thickness, and optical thickness. The sample cloud was located at a height of 1 km with an optical thickness of 2 and a geometrical thickness of 500 m. The liquid water content was represented with 1, 2, 5, or 10 layers, and the corresponding vertical distribution of liquid water varied from a constant value (1 layer) to an approximately adiabatic profile (10 layers). The fourth scenario (Fig. 5d) was a combination of the first two cases comparing a low (height of 500 m) geometrically thick (1 km) cloud with a high (height of 7 km) geometrically thin (200 m) cloud, both having the same optical thickness. As described above, all calculations were performed for a constant cloud optical thickness of 2. For comparison, all calculations were repeated with a constant cloud optical thickness of 20.

Fig. 5.

Schematics of the sensitivity studies on the influence of the vertical distribution of cloud properties on surface shortwave irradiance. (a) Influence of the vertical location of a given cloud—the same (1D) cloud is placed at five different altitudes in the atmosphere, (b) influence of the geometrical thickness of the cloud for a given optical thickness, (c) influence of the vertical profile in the cloud, and (d) combined influence of vertical position and geometrical thickness.

Fig. 5.

Schematics of the sensitivity studies on the influence of the vertical distribution of cloud properties on surface shortwave irradiance. (a) Influence of the vertical location of a given cloud—the same (1D) cloud is placed at five different altitudes in the atmosphere, (b) influence of the geometrical thickness of the cloud for a given optical thickness, (c) influence of the vertical profile in the cloud, and (d) combined influence of vertical position and geometrical thickness.

The results of the calculations, with a solar zenith angle of 0°, are summarized in Table 1, which shows the maximum difference among the surface irradiances calculated within each of the four scenarios, expressed as a percentage of the mean. The differences among the calculated profiles are below 5% for scenarios 1 and 2, below 1% for scenario 4, and below 0.1% for scenario 3. These results imply that for a given optical thickness the vertical distribution of liquid water content has only little influence on surface irradiance in comparison with the variability of surface irradiance under cloudy conditions. Therefore a lookup table approach is accurate enough for our application, which further reduces the computational time. Because the radiative calculation is done offline, we can even afford more accurate solutions than in an online calculation, such as a discrete ordinate model and a detailed correlated-k distribution, both of which would be too expensive in an online calculation.

Table 1.

Results of the four scenarios of the sensitivity studies on the influence of the vertical distribution of cloud properties on surface shortwave irradiance (illustrated in Fig. 5) for two optical thicknesses, τ = 2 and τ = 20. The maximum difference between the surface irradiances calculated within each of the four scenarios, expressed as a percentage of the mean, is specified.

Results of the four scenarios of the sensitivity studies on the influence of the vertical distribution of cloud properties on surface shortwave irradiance (illustrated in Fig. 5) for two optical thicknesses, τ = 2 and τ = 20. The maximum difference between the surface irradiances calculated within each of the four scenarios, expressed as a percentage of the mean, is specified.
Results of the four scenarios of the sensitivity studies on the influence of the vertical distribution of cloud properties on surface shortwave irradiance (illustrated in Fig. 5) for two optical thicknesses, τ = 2 and τ = 20. The maximum difference between the surface irradiances calculated within each of the four scenarios, expressed as a percentage of the mean, is specified.

For the calculation of the lookup tables we used the libRadtran radiative transfer package (Mayer and Kylling 2005). In particular, the discrete ordinate radiative transfer (DISORT) solver by Stamnes et al. (1988) was used to calculate surface solar irradiance as a function of the integrated cloud optical thickness and solar zenith angle using the k distribution of Kato et al. (1999). The cloud has its base at 1000 m, a geometrical thickness as a function of optical thickness, and an adiabatic liquid water content. The effective droplet radius is a function of liquid water content. For the conversion from cloud microphysical to optical properties, the parameterization of Hu and Stamnes (1993) was used. Profiles of pressure, temperature, water vapor, and so on were taken from the midlatitude summer profile by Anderson et al. (1986). Because water vapor absorption affects surface irradiance to some degree, the lookup table could of course be quickly recalculated for any given profile. The spectral albedo of grass (Bowker et al. 1985) was used as surface albedo.

For application in the CRM, the “tilted optical thickness” τ has to be calculated for each surface grid cell from which the corresponding surface irradiance can then immediately be obtained from the lookup table. Optical thickness can be computed from liquid water content from the CRM by computing the effective droplet radius and using a parameterization to convert liquid water content and effective droplet radius to optical thickness like in Hu and Stamnes (1993). To compute the effective radius reff we used the approach in Martin et al. (1994), who computed the effective droplet radius in micrometers as a function of liquid water content wL (g m−3):

 
formula

with density ρ of water (g m−3), cloud droplet number concentration N [which we assumed to be constant (200 × 106 m−3)], and the constant kM, which assumes values of 0.67 for continental air masses and 0.8 for maritime air masses, respectively. Here the average of both values, kM = 0.735, was used. Liquid water content and effective droplet radius were converted to an extinction coefficient βext (m−1) at 550 nm with the parameterization of Hu and Stamnes (1993) (the lookup table tabulates irradiance as a function of visible optical thickness at 550 nm, but the spectral variability of the optical properties is, of course, taken into account in the calculation of integrated solar irradiance):

 
formula

with βext in inverse kilometers, reff in micrometers, wL in grams per meter cubed, and the values for a, b, and c as given in Table 2 

Table 2.

Parameters for the calculation of the extinction coefficient, following the parameterization of Hu and Stamnes (1993).

Parameters for the calculation of the extinction coefficient, following the parameterization of Hu and Stamnes (1993).
Parameters for the calculation of the extinction coefficient, following the parameterization of Hu and Stamnes (1993).

d. Parallelization

The large-eddy simulation model EULAG is parallelized and usually runs on several processors. For that purpose, the model domain is divided horizontally into subdomains, each containing a fraction of the grid points in the horizontal direction. Depending on solar zenith angle, tilted rays in the TICA scheme may cross more than one subdomain. In our implementation, we use one ray for each individual surface grid cell, centered in the grid cell, as representative for the tilted column (in principle, one should average over many rays starting from different locations in the grid cell, but as the validation later in this manuscript shows, one ray is enough for our application). Thus, we need to accumulate the optical thickness for NxNy rays, where Nx and Ny are the numbers of grid points in the horizontal x and y directions, respectively. Integration of the optical thickness along each individual ray would require the information of the cloud properties from all the subdomains that it penetrates. This would be very inefficient in a parallelized code because it requires permanent communication between processors. To circumvent this problem, we rather calculate the contribution of each subdomain to each surface grid cell (i.e., the optical thickness of each surface grid cell accumulated in the subdomain), instead of following each ray through the complete grid. The contributions of all subdomains are added in a final step, and thus communication between processors is required only once for each model time step and only cloudy grid boxes have to be considered.

EULAG defines properties such as liquid water content at grid points, whereas we require grid boxes to accumulate the optical thickness. For that purpose, we define grid boxes around each model grid point within which we assume constant liquid water content. The box (i, j, k) ranges from zk to zk+1, from xi − Δx/2 to xi + Δx/2, and from yj − Δy/2 to yj + Δy/2, where Δx and Δy are the grid spacings in the x and y directions, respectively. The number of rays and the length of the ray penetrating each box depend on sun position. With overhead sun (i.e., solar zenith angle equal to zero), exactly one ray runs through each box. The length of this ray within the box is equal to the vertical extension Δz of the box. With increasing solar zenith angle, the number of rays running through each box increases. Provided that all boxes have the same horizontal extensions Δx and Δy, exactly one ray enters each box at its bottom (see Fig. 6). This ray is considered first: Initially, the index of the surface grid cell that this ray finally hits has to be calculated. The base area of the box (i, j, k) is projected to the base area of the model domain according to the solar zenith θ and azimuth angle ϕ (ϕ is equal to 0° when the sun is in the south). This is shown in Fig. 6 (cross section) and Fig. 7 (top view). The projection of the point (xi − Δx/2, yj − Δy/2, zk) results in the following coordinates:

 
formula
 
formula
Fig. 6.

Determination of the index of the ray entering the current box through the bottom (cross section): the bottom of the box (i, j, k) is projected to the bottom of the model domain along the dashed lines. Rays from the center of the surface pixel in the direction of the sun are indicated with solid lines.

Fig. 6.

Determination of the index of the ray entering the current box through the bottom (cross section): the bottom of the box (i, j, k) is projected to the bottom of the model domain along the dashed lines. Rays from the center of the surface pixel in the direction of the sun are indicated with solid lines.

Fig. 7.

Determination of the index of the ray entering the current box through the bottom (top view): the bottom of the box (i, j, k) is projected to the bottom of the model domain along the dashed lines; the shaded area marks the projected box bottom. A ray reaching from the center of the respective surface pixel in the direction of the sun is indicated with a solid line.

Fig. 7.

Determination of the index of the ray entering the current box through the bottom (top view): the bottom of the box (i, j, k) is projected to the bottom of the model domain along the dashed lines; the shaded area marks the projected box bottom. A ray reaching from the center of the respective surface pixel in the direction of the sun is indicated with a solid line.

Here θ < 90° is required. In the case of the projection being outside the base area of the model domain, it will be shifted into the model domain according to periodic boundary conditions used in the CRM. By comparing the projections of the four corners (xproj,i, yproj,j), (xproj,i+1, yproj,j), (xproj,i+1, yproj,j+1), and (xproj,i, yproj,j+1) with the coordinates of the base pixels, the index (is, js) of the considered ray is obtained. The projection of the current model box covers exactly one central point of a pixel at the base of the model domain.

Using the index (is, js), the interception point (xs, ys) of this ray with the base of the current model box (i, j, k) can be calculated using

 
formula
 
formula

For the calculation of the length of the ray in the current box, the projections lx, ly, and lz of the ray up to the interception point with the next y–z, xz, or xy plane are needed (see Fig. 8). These are calculated for a ray in the positive x and y directions following

 
formula
 
formula
 
formula

For a ray proceeding in the negative x and y directions, lx, ly, and lz are calculated analogously. The length of the ray in the box (i, j, k) can then be obtained using

 
formula

In the case of ϕ equal to 0° or 180°, lx is not defined, because there is no interception with a y–z plane if the ray is running parallel to the x axis. In a similar way, ly is not defined for ϕ equal to 90° or 270°.

Fig. 8.

Calculation of the length of a ray through the box (i, j, k) (top view).

Fig. 8.

Calculation of the length of a ray through the box (i, j, k) (top view).

The contribution of the box (i, j, k) to the optical thickness of the ray is finally calculated using the length lB and the extinction coefficient βext(i, j, k), determined according to Eq. (2):

 
formula

Following the calculation of the optical thickness, the point (xe, ye, ze) at which the ray leaves the box is calculated using

 
formula
 
formula
 
formula

If this point is at the top of the box, no further ray is intercepting this box. In the case of the point (xe, ye, ze) being on one of the sidewalls of the grid box, another ray enters this box at the same height at the opposite wall. This ray will be handled according to Eqs. (7)(14). After calculating the contributions of all cloudy grid boxes within a subdomain, the optical thickness of each surface grid cell accumulated in the subdomain is obtained. In a final step, the contributions of all subdomains are added and the surface irradiance is obtained as a function of the optical thickness using the lookup table.

The parallelization is realized in a way that reduces computational demand considerably. For the simulation of a broken cloud field, the calculation time of the model with the new algorithm is increased by only 3% relative to the original version of the model, which did not include any radiation calculation.

3. Validation

To validate the presented method, the algorithm has been applied to different cloud fields and for different solar zenith and azimuth angles. The results were compared with fully 3D radiative transfer calculations by the Monte Carlo model MYSTIC (Mayer 2000a, b). MYSTIC has been extensively validated in the Intercomparison of 3D Radiation Codes (I3RC; Cahalan et al. 2005), where the irradiance computed by MYSTIC usually agreed to better than 1% with other 3D models.

An example is presented in Figs. 9 and 10. The numerically simulated cloud field used for the comparison is represented in Fig. 9 by its liquid water path. The right panel shows an inclined view of the cloud field. The cloud field with a size of 6.4 km × 6.4 km consists of cumulus mediocris with a cloud cover of 0.3, the cloud base is at 600 m, and the cloud top is at 1600 m. The lines indicate the cross sections selected for comparison: they were chosen to represent different areas with large, middle, and small shaded areas. Figure 10 (top left) shows the nadir radiance distribution at the top of the atmosphere simulated with MYSTIC, to illustrate how one would see this cloud scene from above (e.g., from an airplane). For this calculation, the sun is in the south with a solar zenith angle of 30°. The comparisons (Fig. 10) between the 3D (solid) and TICA with (dashed) and without (dotted) smoothing of the diffuse irradiance show that the locations of cloud shadows are perfectly represented by the TICA calculations. The smoothing of the diffuse irradiance yields a considerable improvement that results in a very good agreement between the expanded TICA and the 3D calculations. The reduction of the irradiance in the shadows is quantitatively reproduced, as is the enhancement of the irradiance outside the shadows (for comparison, the TICA without smoothing is equal to the corresponding cloudless sky value outside the shadows). An overview of the results of the validation is given in Table 3, which shows the results of the new TICA scheme (with smoothing of the diffuse irradiation) as well as the results of the standard TICA (without smoothing) and of the ICA (calculated with the lookup tables). The overestimation or underestimation of the mean surface irradiance relative to 3D calculations and the error, given by σ(ETICAE3D)/mean(E3D) with the standard deviation σ, are shown for solar zenith angles 30° and 60°. The smoothing of the diffuse irradiance reduces the error of the TICA considerably. The error of the ICA is 2 times that of the standard TICA.

Fig. 9.

(left) Vertically integrated liquid water path of a simulated cumulus cloud field (the vertical lines marked with an A, B, and C indicate the location of the calculations shown later in Fig. 10). (right) Inclined view on this cloud field (gray tone indicates liquid water content, and the labeling of the axes corresponds to the model grid).

Fig. 9.

(left) Vertically integrated liquid water path of a simulated cumulus cloud field (the vertical lines marked with an A, B, and C indicate the location of the calculations shown later in Fig. 10). (right) Inclined view on this cloud field (gray tone indicates liquid water content, and the labeling of the axes corresponds to the model grid).

Fig. 10.

(top left) Reflectivity of a simulated cumulus cloud field. Also shown is a comparison of surface irradiance between 3D (solid line) and TICA with (dashed line) and without (dotted line) smoothing of diffuse irradiance along the vertical lines shown in the left panel of Fig. 9: (top right) line A, (bottom left) line B, and (bottom right) line C.

Fig. 10.

(top left) Reflectivity of a simulated cumulus cloud field. Also shown is a comparison of surface irradiance between 3D (solid line) and TICA with (dashed line) and without (dotted line) smoothing of diffuse irradiance along the vertical lines shown in the left panel of Fig. 9: (top right) line A, (bottom left) line B, and (bottom right) line C.

Table 3.

Comparison between 3D radiative transfer calculations and the new TICA algorithm and the standard TICA and ICA for a cumulus mediocris (stratocumulus) cloud field.

Comparison between 3D radiative transfer calculations and the new TICA algorithm and the standard TICA and ICA for a cumulus mediocris (stratocumulus) cloud field.
Comparison between 3D radiative transfer calculations and the new TICA algorithm and the standard TICA and ICA for a cumulus mediocris (stratocumulus) cloud field.

The same analysis was repeated for overcast clouds. The overcast stratocumulus cloud field (not shown) had mean optical and geometrical thicknesses of 26 and 400 m, respectively. The results are as shown in Table 3. For this overcast cloud field, the errors for both the ICA and the TICA are much smaller. As expected, the TICA is mainly required for broken cloud fields with shadows; for more homogeneous, overcast cases, the ICA produces reasonable results. However, even in the overcast case, the smoothing of the diffuse irradiance reduces the uncertainty considerably, illustrating that the ICA produces too-sharp contrasts.

4. Summary and conclusions

A fast method has been developed to include accurate 3D radiative transfer calculations in an LES, based on the tilted independent column approximation. By smoothing the diffuse irradiance, we obtained very high accuracy relative to real 3D calculations. The method was validated by comparison with a detailed 3D radiation code and was implemented into the large-eddy simulation model EULAG.

In contrast to the widely employed independent column approximation, the TICA correctly reproduces not only mean values but also the spatial distribution of the integrated shortwave surface flux. Whereas the ICA should only be applied for overcast cloud situations or solar zenith angles close to zero, the TICA can be used for any cloud field and arbitrary solar zenith angles. The possibility to parallelize this algorithm enhances its applicability in CRMs. In our example implementation, the computational time in EULAG was increased by only 3% relative to the reference run without radiation. In a follow-on paper, we will present a detailed study of the effect of cloud shadows on cloud formation and evolution through the feedback of reduced solar irradiance on surface heating.

REFERENCES

REFERENCES
Anderson
,
G.
,
S.
Clough
,
F.
Kneizys
,
J.
Chetwynd
, and
E.
Shettle
,
1986
:
AFGL atmospheric constituent profiles (0–120 km). Tech. Rep. AFGL-TR-86–0110, Hanscom AFB, MA, 46 pp
.
Bowker
,
D. E.
,
R. E.
Davis
,
D. L.
Myrick
,
K.
Stacy
, and
W. T.
Jones
,
1985
:
Spectral reflectances of natural targets for use in remote sensing studies. NASA Ref. Publ. 1139, 188 pp
.
Cahalan
,
R. F.
,
W.
Ridgway
, and
W.
Wiscombe
,
1994
:
Independent pixel and Monte Carlo estimates of stratocumulus albedo.
J. Atmos. Sci.
,
51
,
3776
3790
.
Cahalan
,
R. F.
, and
Coauthors
,
2005
:
The I3RC: Bringing together the most advanced radiative transfer tools for cloudy atmospheres.
Bull. Amer. Meteor. Soc.
,
86
,
1275
1293
.
Grabowski
,
W. W.
, and
P. K.
Smolarkiewicz
,
2002
:
A multiscale anelastic model for meteorological research.
Mon. Wea. Rev.
,
130
,
939
956
.
Hu
,
Y. X.
, and
K.
Stamnes
,
1993
:
An accurate parameterization of the radiative properties of water clouds suitable for use in climate models.
J. Climate
,
6
,
728
742
.
Kato
,
S.
,
T. P.
Ackerman
,
J. H.
Mather
, and
E. E.
Clothiaux
,
1999
:
The k-distribution method and correlated-k approximation for a shortwave radiative transfer model.
J. Quant. Spectrosc. Radiat. Transfer
,
62
,
109
121
.
Markowski
,
P. M.
, and
J. Y.
Harrington
,
2005
:
A simulation of a supercell thunderstorm with emulated radiative cooling beneath the anvil.
J. Atmos. Sci.
,
62
,
2607
2617
.
Marshak
,
A.
,
A.
Davis
,
R.
Cahalan
, and
W.
Wiscombe
,
1998
:
Nonlocal independent pixel approximation: Direct and inverse problems.
IEEE Trans. Geosci. Remote Sens.
,
36
,
192
204
.
Martin
,
G. M.
,
D. W.
Johnson
, and
A.
Spice
,
1994
:
The measurement and parameterization of effective radius of droplets in warm stratocumulus clouds.
J. Atmos. Sci.
,
51
,
1823
1842
.
Mayer
,
B.
,
2000a
:
I3RC phase 1 results from the MYSTIC Monte Carlo model.
Intercomparison of Three-Dimensional Radiation Codes: Abstracts of the First and Second International Workshops, R. F. Cahalan and R. Davies, Eds., Institute of Atmospheric Physics, University of Arizona, 49–54
.
Mayer
,
B.
,
2000b
:
I3RC phase 2 results from the MYSTIC Monte Carlo model.
Intercomparison of Three-Dimensional Radiation Codes: Abstracts of the First and Second International Workshops, R. F. Cahalan and R. Davies, Eds., Institute of Atmospheric Physics, University of Arizona, 107–108
.
Mayer
,
B.
, and
A.
Kylling
,
2005
:
Technical note: The libRadtran software package for radiative transfer calculations: Description and examples of use.
Atmos. Chem. Phys.
,
5
,
1855
1877
.
Mims
,
F.
, and
J.
Frederick
,
1994
:
Cumulus clouds and UV-B.
Nature
,
371
,
291
.
O’Hirok
,
W.
, and
C.
Gautier
,
1998
:
A three-dimensional radiative transfer model to investigate the solar radiation within a cloudy atmosphere. Part I: Spatial effects.
J. Atmos. Sci.
,
55
,
2162
2179
.
Schumann
,
U.
,
A.
Dörnbrack
, and
B.
Mayer
,
2002
:
Cloud-shadow effects on the structure of the convective boundary layer.
Meteor. Z.
,
11
,
285
294
.
Smolarkiewicz
,
P. K.
, and
L. G.
Margolin
,
1997
:
On forward-in-time differencing for fluids: An Eulerian/semi-Lagrangian non-hydrostatic model for stratified flows.
Atmos.–Ocean
,
35
,
127
152
.
Stamnes
,
K.
,
S.
Tsay
,
W.
Wiscombe
, and
K.
Jayaweera
,
1988
:
A numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media.
Appl. Opt.
,
27
,
2502
2509
.
Várnai
,
T.
, and
R.
Davies
,
1999
:
Effects of cloud heterogeneities on shortwave radiation: Comparison of cloud-top variability and internal heterogeneity.
J. Atmos. Sci.
,
56
,
4206
4224
.
Weinman
,
J.
, and
Harshvardhan
,
1982
:
Solar reflection from a regular array of horizontally finite clouds.
Appl. Opt.
,
21
,
2940
2944
.

Footnotes

* Current affiliation: School of Earth Sciences, The University of Melbourne, Melbourne, Australia.

Corresponding author address: Kathrin Wapler, Deutscher Wetterdienst/German Weather Service, Frankfurterstr. 135, 63067 Offenbach, Germany. Email: kathrin.wapler@dwd.de