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%.
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 m−2. 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.
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).
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.
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.
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.
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.
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.
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):
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):
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
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:
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
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, x–z, or x–y plane are needed (see Fig. 8). These are calculated for a ray in the positive x and y directions following
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
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°.
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):
Following the calculation of the optical thickness, the point (xe, ye, ze) at which the ray leaves the box is calculated using
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.
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 σ(ETICA − E3D)/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.
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.
* 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: firstname.lastname@example.org