Quantifying the bias of radiative heating rates in NWP models for shallow cumulus clouds

The interaction between radiation and clouds represents a source of uncertainty in numerical weather prediction (NWP) due to both intrinsic problems of one-dimensional radiation schemes and poor representation of clouds. The underlying question addressed in this study is how large is the NWP radiative bias for shallow cumulus clouds and how does it scale with various input parameters of radiation schemes, such as solar zenith angle, surface albedo, cloud cover and liquid water path. A set of radiative transfer calculations was carried out for a realistically evolving shallow cumulus cloud field stemming 5 from a large-eddy simulation (LES). The benchmark experiments were performed on the highly-resolved LES cloud scenes using a three-dimensional Monte Carlo radiation model. An absence of middle and high cloud is assumed above the shallow cumulus cloud layer. In order to imitate poor representation of shallow cumulus in NWP models, cloud optical properties were horizontally averaged over the cloudy part of the boxes with dimensions comparable to NWP horizontal grid spacing (several km) and the common δ-Eddington two-stream method with maximum-random overlap assumption for partial cloudiness was 10 applied (denoted as "1-D" experiment). The bias of the 1-D experiment relative to the benchmark was investigated in the solar and thermal part of the spectrum, examining the vertical profile of heating rate within the cloud layer and net surface flux. It is found that during daytime and nighttime, the destabilization of the cloud layer in the benchmark experiment is artifically enhanced by an overestimation of the cooling at cloud top and an overestimation of the warming at cloud bottom in the 1-D experiment (bias of about −15 K day−1 is observed locally for stratocumulus scenarios). This destabilization, driven by the 15 thermal radiation, is maximized during nighttime, since during daytime the solar radiation has a stabilizing tendency. The daytime bias at the surface is governed by the solar fluxes, where the 1-D solar net flux overestimates (underestimates) the corresponding benchmark at low (high) sun. The overestimation at low sun (bias up to 80 % over land and ocean) is largest at intermediate cloud cover, while underestimation at high sun (bias up to −40 % over land and ocean) is peaked at larger cloud cover (80 % and beyond). At nighttime, the 1-D experiment overestimates the amount of benchmark surface cooling 20 with the maximal bias of about 50 % peaked at intermediate cloud cover. Moreover, an additional experiment was carried out by running the Monte Carlo radiation model in the independent column mode on cloud scenes preserving their LES structure (denoted as "ICA" experiment). The ICA is predominantly more accurate than the 1-D experiment (with respect to the same benchmark). This highlights the importance of an improved representation of clouds even at the resolution of today’s regional (limited-area) numerical models, which needs to be considered if NWP radiative biases are to be efficiently reduced. All in all, 25 1 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2018-1247 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 18 January 2019 c © Author(s) 2019. CC BY 4.0 License.


Introduction
Despite great progress in the availability of computing power over the recent years, exact three-dimensional (3-D) radiative transfer solvers remain computationally too expensive to be routinely used in numerical weather prediction (NWP). Therefore 5 approximate one-dimensional (1-D) solvers are applied operationally in NWP models, introducing uncertainty to weather forecasts. The underlying assumption of 1-D solvers is the Independent Column Approximation (ICA), where the radiative transfer problem is solved in each model column separately, thus entirely neglecting horizontal photon transport between individual columns across the model domain. As a standard technique for the solution of the radiative transfer equation within an ICA column computationally efficient analytical two-stream methods (TSMs) are widely employed, having a laudable 10 tradition starting more than a century ago with Schuster (1905) and Schwarzschild (1906).
A radiation parameterization scheme within a host NWP model has a twofold task: Apart from the calculation of radiative heating rate distribution within the atmosphere [K day −1 ] for the diabatic heating term in the prognostic temperature equation, it needs to supply the net surface radiative flux [W m −2 ] for the proper evaluation of the surface energy budget. While in clear-sky conditions 1-D solvers generally provide good estimates for atmospheric heating rates and surface fluxes, the error 15 related to neglected 3-D effects considerably increases in the presence of clouds (e.g., Schmetz, 1984;O'Hirok and Gautier, 1998a, 1998b, 2004Tompkins, 2003, 2005;Tompkins and Di Giuseppe, 2007;Wissmeier et al., 2013).
Although the uncertainty related to the inaccurate treatment of radiation-cloud interaction in numerical models depends, inter alia, on cloud type, the present study is restricted to shallow cumulus clouds. The importance of shallow convection for the redistribution of atmospheric heat and moisture is well-acknowledged (e.g., Tiedtke, 1989;Albrecht 20 et al., 1995;Zhao and Austin, 2005). Small-scale fluctuations in microphysical, dynamical and thermodynamical parameters observed in boundary layers containing cumuli (Barker et al., 1985) indicate the complexity of a competent coupling of cumulus cloud fields with a full 3-D radiative field.
In the last decades large-eddy simulation (LES) and cloud-resolving (CRM) models have become an important tool in boundary layer research (Neggers et al., 2003). If an accurate 3-D radiation (e.g., Monte Carlo) calculation is performed on 25 such 3-D highly-resolved cumulus cloud ("benchmark experiment" or "truth"), the following is observed: In the solar spectral range the largest heating is at the illuminated cloud side and the shadow of the cloud at the ground is shifted according to solar zenith angle (Wapler and Mayer, 2008;Wissmeier et al., 2013;Mayer, 2015, 2016). In the thermal spectral range there is a strong cooling of cloud top and cloud sides and modest warming of cloud bottom (Kablick III et al., 2011;Mayer, 2014, 2016). In the associated ICA calculation, the 3-D effects related to cloud sides are misrepresented, the main 30 shortcomings are as follows: In the solar spectral range the heating is always at cloud top and the shadow at the surface lies directly underneath the cloud, which is fundamentally wrong unless the sun is in the zenith Mayer, 2015, 2016).
In the thermal spectral range, the ICA approximation only captures cloud top cooling and cloud base warming, but entirely neglects cooling of cloud sides (Kablick III et al., 2011;Mayer, 2014, 2016).
When 3-D radiation parameterizations are coupled to a LES or CRM model, the above-mentioned shortcomings of ICA on the cumulus cloud evolution can be studied. Klinger et al. (2017) showed that interactive 3-D thermal radiation affects cloud circulation by enhancing cloud-core updrafts and surrounding subsiding shells. In addition, it alters the organization of 5 clouds (i.e., convective self-aggregation). In a recent study by Jakub and Mayer (2017) it was shown that interactive 3-D solar radiative transfer may cause formation of cloud streets similar to the known roll clouds caused by wind shear, whereas the ICA approximation produces randomly positioned clouds.
Before these effects can be directly simulated within NWP, we are presumably at least a decade away from the desired resolution. Today's regional NWP models with horizontal grid spacing on the order of few km (e.g., operational model of Ger-10 man Weather Service in its convection-permitting configuration COSMO-DE with horizontal grid spacing of 2.8 km) mostly resolve deep convection and have a parameterization for shallow convection. Depending on cloud parameterization scheme subgrid-scale cloudiness (cloud fraction) within a model grid box is usually diagnosed from the grid-scale relative humidity.
Shallow cumulus clouds in state-of-the-art NWP are thus represented as horizontally homogeneous layers of partial cloudiness, entirely missing their 3-D geometrical structure and small-scale variability of optical properties. Further, an assumption is re- 15 quired of how partial cloudiness is distributed in the vertical direction, which is another deficiency of NWP radiation schemes (e.g., Barker et al., 2003;Barker, 2008;Wu and Liang, 2005;Kablick III et al., 2011). The widely employed assumption is the maximum-random overlap assumption (Geleyn and Hollingsworth, 1979), which, implemented in the two-stream framework, gives the so called two-stream method with maximum-random overlap assumption for partial cloudiness (Ritter and Geleyn, 1992), commonly used as NWP radiation solver.

20
To summarize, the interaction between radiation and shallow cumulus clouds represents a source of uncertainty in NWP due to both intrinsic problems of 1-D radiation schemes and poor representation of subgrid-scale clouds. The underlying question of the present study is as follows: How large is the bias of NWP radiative heating rates on shallow cumulus clouds and how does it scale with various input parameters of radiation schemes, such as solar zenith angle (SZA), surface albedo (A), total cloud cover (CC) and cloud liquid water path (LWP)? 25 This paper is organized as follows: The experimental design is outlined in section 2, whereas in the subsequent section 3 the results are presented. Summary and concluding remarks are given in section 4.

Theory and method
First, in section 2.1, the Monte Carlo radiation model, which is used for the 3-D benchmark and ICA calculations, is introduced.
Then, in section 2.2, we describe the δ-Eddington two-stream method with maximum-random overlap assumption for partial 30 cloudiness as employed in this study. In the subsequent sections 2.3 and 2.4 the cloud field data set and the related experimental strategy are presented. Finally, the general setup of radiative transfer computations is summarized in section 2.5.

Reference model -MYSTIC
The benchmark experiments are performed with the 3-D radiative transfer model MYSTIC, the Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres (Mayer, 2009), which is part of the libRadtran software package (www.libradtran.org; Mayer and Kylling, 2005;Emde et al., 2016) and can be run in independent column mode as well.
MYSTIC participated in both phases of the international Intercomparison of 3-D Radiation Codes (I3RC, Cahalan et al.,5 2005), where it proved its ability to accurately compute radiative transfer in versatile cloud scenarios.
2.2 δ-Eddington two-stream method with maximum-random overlap assumption for partial cloudiness We begin by introducing the classic δ-Eddington two-stream method (Joseph et al., 1976) suitable for a horizontally homogeneous model atmosphere and then explain the extension of this method which accounts for partial cloudiness.
The common feature of two-stream methods is the division of the radiation field into direct (unscattered) solar beam (S) 10 and two streams of diffuse radiation − the downward (E ↓ ) and upward (E ↑ ) component. For most applications, δ-TSMs, in which a part of the scattered radiation is retained in the direct beam to approximate the strong forward-scattering peak of cloud droplets and aerosol particles, have been found to be more accurate than two-stream methods without δ-scaling (Räisänen, 2002). The fractional scattering into the forward peak taken to be the square of the phase function asymmetry parameter is what distinguishes the widely used δ-Eddington approximation from others of similar nature. 15 For the calculations in a vertically inhomogeneous atmosphere, the atmosphere is discretized into a number of homogeneous layers, each characterized by its optical properties (optical thickness, single scattering albedo, asymmetry parameter). Consider first a single layer (j) located between levels (i-1) and (i) 1 . A system of linear equations determining the fluxes that emanate from this layer as a function of fluxes that enter the layer can be written in matrix form as: a 11 a 12 a 13 a 12 a 11 a 23 0 0 a 33 (1)

20
The linear coefficients a kl in Eq. (1), referred to as Eddington coefficients, are functions of optical properties of layer (j). They have the following physical meaning: a 11 and a 12 represent the transmission and reflection coefficient for diffuse radiation, respectively. Further, a 13 and a 23 represent the reflection and transmission coefficient for the primary scattered solar radiation, respectively, while a 33 denotes the transmission coefficient for the direct solar radiation. For the details of their definition, as well as for the inclusion of thermal radiation, the reader is referred to Zdunkowski (1980). 25 In a consecutive step individual layers are concatenated by imposing flux continuity at each level. Taking appropriate boundary conditions at the top of the atmosphere (TOA) and at the ground into account, the equation system is solved analytically by means of standard numerical procedures (e.g., Zdunkowski, 1980;Stephens and Webster, 1981;Ritter and Geleyn, 1992; 1 We follow the convention of i, j increasing downward from the top of the atmosphere, where i = 0, j = 1. Index i is used for level variables, whereas index j is used for layer variables. The N vertical layers, which are enumerated from 1 to N , are enclosed by (N +1) vertical levels, which are enumerated from 0 to N . Stephens et al., 2001). After radiative fluxes thrughout the atmosphere have been computed, the calculation of heating rates is straightforward. The heating rate [K day −1 ] of an individual layer is given by: where ρ represents air density, c p represents the specific heat capacity of air at constant pressure, ∆z represents the vertical thickness of the layer and ∆E net represents the radiative flux absorbed in the layer [W m −2 ], defined as the difference between 5 the fluxes entering the layer and those leaving the layer.
Consider now a partially cloudy layer (Fig. 1, left panel), which is characterized by two sets of optical properties and corresponding Eddington coefficients − one for the cloudy region (superscript c) and the other for the cloud-free region (superscript f ). In order to apply the maximum-random overlap assumption (Geleyn and Hollingsworth, 1979) the cloudy and cloud-free fluxes need to be treated separately. Total radiative flux at a given level is thus the sum of the cloudy and cloud-free 10 component, e.g., and analogously for both diffuse components. The equation system (1) is replaced by: 15 so that the fluxes emanating from the cloudy and cloud-free region depend on a linear combination of both cloudy and cloudfree incoming fluxes. Overlap coefficients p 1 , p 2 , p 3 and p 4 refer to layer (j) and describe the division of incoming fluxes between the cloudy and cloud-free region in accordance with the maximum-random overlap assumption, where adjacent cloudy layers are overlapped maximally, and cloudy layers separated by at least one cloud-free layer are overlapped randomly. For layer (j) they have the following form: p 3 (j) = min C(j), C(j − 1) p 4 (j) = min C(j), C(j + 1) with C representing layer cloud fraction. Figure 1 (right and middle panel) illustrates both possible geometries that need to be considered in order to determine the coefficients p 1 and p 3 related to the division of downward fluxes (the division of upward fluxes is managed via p 2 and p 4 in a similar fashion). The method has been successfully implemented in the radiative transfer 5 package libRadtran for the purpose of this study.

Cloud field data set
Input for offline radiation calculations is a set of shallow cumulus cloud fields, simulated with the University of California, Los (1) in section 3.1.3). For our analysis we choose a set of ten cloud scenes (depicted in Fig. 2) from the initial 8 hours of simulation in a way that total cloud cover of the scene is varied between ∼10 % and ∼100 % approximately with a step of 10 %. Thus the set comprises examples of broken cumulus as well as more uniform 15 stratocumulus clouds. These cloud scenes have highly variable optical thickness, with a maximum vertically integrated optical thickness of ∼20, ∼80 and ∼230 corresponding to scenes with total cloud cover of ∼10 %, ∼50 % and ∼100 %, respectively.

Strategy
To assess the bias of NWP radiative quantities, an NWP-type experiment together with a benchmark is required. The benchmark calculation using MYSTIC was performed on a cloud field preserving its LES resolution, with the result horizontally averaged over the domain (hereafter abbreviated to "3-D" experiment). In order to mimic poor representation of shallow cumulus in NWP models and thus create a proper NWP-type experiment, the information content of the cloud field needs to be reduced. Therefore 5 LWC and R e were horizontally averaged over cloudy part of the boxes with dimensions comparable to NWP horizontal grid spacing (3.2 km). In this way, four NWP-sized boxes which generally contain partial cloudiness were created in each vertical layer ( Fig. 3, middle panel) and the δ-Eddington two-stream method with maximum-random overlap assumption was called four times per cloud scene. The resulting radiative quantities were again horizontally averaged over the domain (hereafter abbreviated to "1-D" experiment). Moreover, it should be noted that the resolution of the cloud field in the 1-D experiment 10 was only degraded in the horizontal plane, whereas the vertical resolution was kept the same as inherited from LES grid (25  conditions in NWP models and serve as input for two-stream method with maximum-random overlap assumption. The third panel from the left shows total cloud cover versus maximal layer cloud fraction for the entire set of ten cumulus cloud scenes. The discrepancy between the two curves is a measure of the erroneousness of the maximum-random overlap assumption. The right panel shows vertical profiles of cloud fraction for the three selected cloud scenes with total cloud cover of 12.0 %, 52.3 % and 98.9 %. m). 2 In summary, the 1-D experiment has multiple error sources, namely, the poor horizontal cloud structure, vertical overlap assumption, as well as the neglected grid-scale and subgrid-scale horizontal photon transport.
Furthermore, we created a third, intermediate experiment by running MYSTIC in independent column mode on a cloud field preserving its LES resolution and again averaging the result horizontally over the domain in the final step (hereafter abbreviated to "ICA" experiment). The ICA experiment is thus the same as the 3-D experiment except that it neglects horizontal photon 5 transport between the LES-grid columns. In this way we are able to isolate and quantify the contribution of neglected horizontal photon transport to the overall error of 1-D experiments (in the hypothetical case when the subgrid-scale cloud structure would be perfectly guessed).
For each of the three experiments we diagnosed radiative heating rate within the cloud layer and net surface flux (i.e., difference between the total downward and upward surface flux). The error is given by the absolute and relative bias and the 10 root mean square error (RMSE): where y denotes the result of either the 1-D or the ICA experiment and x denotes the result of the 3-D experiment. Three-dimensional distributions of LWC and R e were converted into optical properties using parameterization of Hu and Stamnes (1993), which uses the Henyey-Greenstein phase function as an approximation of the real Mie phase function. The background profiles of atmospheric pressure, temperature, density and trace gases (water vapour, O 3 , CO 2 ) were taken from 10 the US standard atmosphere (Anderson et al., 1986) and are horizontally homogeneous across the domain. Parameterization of absorption and scattering properties of the atmosphere in the solar part of the spectrum follows the correlated-k distribution of Kato et al. (1999). Parameterization of molecular absorption in the thermal spectral range was adopted from Fu and Liou (1992). Solar zenith angle was varied from 0 • (overhead sun) to 80 • with a step of 10 • , while solar azimuth angle was held fixed at 0 • (i.e., sun illuminating the cloud scenes from the south). Further, we varied surface shortwave albedo by applying 15 constant values of 0.25 and 0.05 as typical high and low values representing land and ocean. In the thermal part of the spectrum surface was assumed to be non-reflective. In the Monte Carlo calculations in the solar spectral range we applied a standard forward photon tracing method, starting a total number of 65536000 photons (1000 photons per LES-grid column) at TOA. In the thermal spectral range the backward photon tracing method of Klinger and Mayer (2014) was used and 1000 photons were required in each LES-grid box of the cloud layer. additionally the dependence on surface albedo, cloud cover and liquid water path is discussed.
3.1 Heating rate in the cloud layer 25 The vertical profile of radiative heating rate influences atmospheric stratification and directly impacts flow dynamics. Figure 4 shows the vertical profile of radiative heating rate in the trio (3-D, 1-D, ICA) of experiments in the solar spectral range for SZA of 0 • , 30 • and 60 • as well as in the thermal spectral range. In order to highlight the effects of clouds on radiative biases we examine only the profiles within the cumulus cloud layer, which is located between 1.025 and 1.875 km height (the cloud-free atmosphere below and above these heights is not shown). In the solar 3-D experiment for overhead sun (top left panel of Fig. 4) there is strong absorption of solar radiation in the cloud layer, resulting in a peak heating rate of about 9 K day −1 reached at approximately 1.5 km height. This is slightly above the height of maximum cloud fraction (1.4 km) due to the fact that cloud liquid water, which is the dominant absorber of solar radiation in this layer, generally increases with height from cloud base towards cloud top (except in the uppermost region of the cloud, where it decreases with height due to entrainment). The 5 maximum heating rate is thus located between the height of maximum cloud fraction and the height of maximum LWC. As the sun descends, the incoming radiation at TOA decreases with the cosine of SZA and so does the solar heating rate in the cloud layer, reaching a maximum value of about 8 K day −1 at SZA of 30 • and about 5 K day −1 at SZA of 60 • in the 3-D experiments. The height where the maximum heating is reached stays approximately the same for all SZAs. In the thermal spectral range the cloud layer is subjected to strong cooling attaining a peak value of about 13 K day −1 again at a height of 10 ∼1.5 km. Below this height, the magnitude of cooling decreases towards the cloud base, where a slight cloud base warming effect is observed.
The difference between the 1-D experiment and the 3-D benchmark (Fig. 4, right panels) is as described in the following: In the solar experiments at high sun (SZA between 0 • and 50 • ) the bias of the 1-D profile shows pronounced vertical gradient within the cloud layer and changes its sign approximately at a height of maximum cloud fraction. In the top part of the cloud 15 layer (i.e., above maximum cloud fraction) the 1-D solar heating rate is too high, while in the bottom part of the cloud layer it is too low, compared to the 3-D heating rate. In the case of low sun (SZA of 60 • and larger) the 1-D solar heating rate is systematically too low compared to its 3-D counterpart throughout the entire cloud layer. The main reason for that is cloud side illumination Mayer, 2015, 2016), which is taken into account in 3-D experiments and completely absent in 1-D calculations.

20
The ICA calculations help to explain these findings: In the case of overhead sun the amount of radiative energy hitting cloud tops is practically the same in both ICA and 3-D experiment, yet in the 3-D experiment some of the photons escape through cloud sides (an effect known as "cloud side escape or loss of photons"; O'Hirok and Gautier, 1998a; Hogan and Shonk, 2013), whereas in the ICA approximation the photons remain trapped within individual atmospheric columns. As a consequence, the ICA heating rate is larger than the 3-D heating rate within the top part of the cloud layer (top row of Fig. 4). This leakage of  In the thermal spectral range, the 1-D experiment overestimates the amount of 3-D cooling in the top part of the cloud layer, whereas in the bottom part of the cloud layer the situation is reversed (a difference of more than 5 K day −1 between 1-D and 3-D is observed locally). The ICA experiment, however, underestimates the magnitude of 3-D thermal cooling throughout 5 the entire vertical extent of the cloud layer (with a maximum difference of less than 1.5 K day −1 ). The latter is as expected, a manifestation of cooling due to horizontal emission of radiation through cloud side areas, which is suppressed in the ICA approximation, but present in the 3-D benchmark Mayer, 2014, 2015). Altogether, this implies that the error arising from the poor representation of cloud horizontal and vertical structure additionally affecting the 1-D experiment acts to reduce the positive difference between ICA and 3-D (and even turning it into a negative one) within the top part of the

Net surface flux
Net surface radiative flux is directly related to surface heating and thereby affects the development of convection. Furthermore, it enters the surface layer parameterization scheme (e.g., soil and vegetation scheme) of a host NWP model and influences various physical processes therein (e.g., hydrological processes, such as melting of snow). Operational 2-m temperature predictions, moreover, are among forecast products of most interest for users, yet they are still subjected to substantial and consistent 25 regional biases in NWP worldwide, partially arising directly from biases of surface radiative fluxes.
Motivated by the desire to understand the causes of the latter, we explore here the net surface flux in the trio (3-D, 1-D, ICA) of experiments on a single cumulus cloud scene with total cloud cover of 52.3 % over land. In the solar spectral range  Again, we aim to untangle this bias dependence of 1-D experiments on SZA, by first explaining purely the effects of neglected horizontal photon transport. Due to the aforementioned loss of photons through cloud sides, the diffuse downward radiation at the surface in 3-D is larger than in ICA (Wapler and Mayer, 2008). This is the main reason why solar net surface flux in the ICA experiment is underestimated relative to 3-D at sun angles between 0 • and 60 • . At SZAs larger than 60 • , however, the so-called "elongated shadow effect" (Wissmeier et al., 2013), which is generally present for all sun positions 5 except for overhead sun, becomes dominant and solar net surface flux in the ICA experiment is overestimated relative to 3-D.
This is essentially the cloud side illumination effect, where the effective total cloud cover (Di Giuseppe and Tompkins, 2003;Tompkins and Di Giuseppe, 2007;Hinkelman et al., 2007) increases with descending sun and hence also the size of the shadow increases with decreasing solar elevation, which is not taken into account in the ICA. This leads to a considerably reduced direct radiation reaching the ground and thus solar net flux in the 3-D experiment when the sun is lower in the sky.

10
Both aforementioned shortcomings of ICA manifest themselves in the 1-D experiment as well. For overhead sun, the apparent reduction of total cloud cover in the 1-D experiment due to the overlap assumption (maximal layer cloud fraction of 37.2 %, which is effectively the total cloud cover in the 1-D experiment, is appreciably lower than the total cloud cover of the 3-D cloud field, that is 52.3 %, see the third panel of Fig. 3) acts to increase the amount of direct radiation reaching the surface (compared to 3-D or ICA case), which in turn reduces the net surface flux bias in the 1-D experiment (compared to that in the 15 ICA experiment). When the sun is from the side, the effective total cloud cover in the 1-D experiment remains 37.2 %, while in the 3-D experiment it is increased well beyond 52.3 %, which further increases the discrepancy in cloud shadow area at the surface in the 1-D and 3-D experiment. In particular, at SZAs larger than 60 • , both the absolute and relative bias of the 1-D experiment are by at least a factor of 2 larger than the corresponding biases of the ICA experiment.
In the thermal spectral range the surface cools by emitting more radiation than it receives with a net

Dependence on surface albedo
It is well known that NWP models can have large temperature errors at coastlines (Hogan and Bozzo, 2015). Due to their high computational cost, the radiation schemes are often applied on a coarser spatial grid (compared to the grid of NWP dynamical core). In regions along coastlines this implies that radiative quantities computed over the ocean are being used at nearby land grid points (where surface temperature and albedo are very different), or vice versa. The alternative practice is averaging input 10 to the radiation scheme onto the coarser grid, which has similar disadvantages.
In the previous sections 3.1 and 3.2 experiments for a solar surface albedo of 0.25 were presented. Here we discuss how the results of these experiments change when the albedo is reduced to a typical oceanic value of 0.05, focusing on the comparison between 1-D and 3-D quantities. As albedo is thus reduced, the solar heating rate within the cloud layer in the benchmark 3-D experiments is reduced, since less radiation is reflected from the surface back towards the cloud. This reduction is largest in 15 the lower part of the cloud layer, whereas at the height where maximum heating is reached (and above), the albedo effect is only marginal. Further, the reduction of solar heating rate with a decreased albedo is largest when the sun is overhead, where a maximum difference of about 0.8 K day −1 is observed between 3-D heating rate profiles in the experiments with A of 0.25 and 0.05, and reduces in significance with descending sun (at SZA of 80 • , this difference is imperceptible, essentially less than 0.05 K day −1 throughout the entire vertical extent of the cloud layer). This implies that the variation of surface albedo has a 20 comparatively large effect on the benchmark heating rate profile at small SZAs and becomes less important with decreasing elevation of the sun.
More relevant for this study, the difference between the 1-D and 3-D heating rate profile in the calculations with A of 0.05 stays practically the same as in those with A of 0.25 (with a deviation being mostly less than 0.1 K day −1 at all SZAs). This means that the relative bias of the 1-D heating rate profile is generally increased as albedo is decreased. This increase of relative 25 bias is smallest when the sun is overhead and gains in significance with descending sun.
Further, at the surface (Fig. 6), solar net flux in the 3-D experiment is generally increased as albedo is decreased from 0.25 to 0.05 (lower surface reflectivity implies more radiation is absorbed in the surface). This increase is largest for overhead sun,  an additional 6 W m −2 . On the whole, the relative bias of the 1-D net surface flux over the ocean stays approximately the same as over land (with a deviation of less than 2 % at all SZAs).
Referring back to the findings of Hogan and Bozzo (2015), our results regarding the sensitivity of the 3-D benchmarks on the variation of surface albedo confirm that along coastlines the radiative quantities should be computed on the regular grid (at least at higher sun elevations). The fact that the absolute bias of cloud-layer heating rate is approximately the same over 5 land and ocean and that the relative bias of net surface flux over land and ocean is approximately the same as well, indicate the possibility of eliminating one parameter (namely the surface albedo) when developing a correction for NWP radiative quantities (after the robustness of the results is proved for diverse cloud scenarios).

Dependence on cloud cover
We examine now the dependence of heating rates and surface fluxes on cloud cover (in addition to their dependence on SZA 10 and albedo) by analysing the entire data set of ten cumulus cloud scenes. We present the benchmark 3-D experiments first and then discuss the bias of 1-D experiments (noting that the ICA experiments are investigated in the next section 3.5, where additionally the dependence on LWP is examined). Thus in the 3-D solar experiments over land, the heating rate within the cloud layer generally becomes larger with increasing CC of the cloud scene. At overhead sun, for example, a peak heating rate of ∼3 K day −1 , ∼9 K day −1 and ∼43 K day −1 in the experiments on cloud scenes with CC of 12.0 %, 52.3 % and 98.9 %, 15 respectively (Fig. 7, top left panel). The height where the peak heating is reached does not vary with SZA and is slightly above the height of maximum cloud fraction of a given cloud scene (the latter differs from scene to scene, since both the vertical extent of the cumulus cloud layer as well as the height of the maximum cloud fraction generally increase during the course of UCLA-LES simulation). Similarly, in the 3-D thermal experiments, the main cloud-radiative effects described in section 3.1 (cloud top cooling, cloud side cooling, cloud base warming) in general become more pronounced as CC is increased. The peak 20 magnitude of cooling, for example, equals to ∼4 K day −1 , ∼13 K day −1 and ∼76 K day −1 in the experiments on cloud scenes with CC of 12.0 %, 52.3 % and 98.9 %, respectively (Fig. 7, bottom left panel). The height where this peak thermal cooling is attained at a given cloud scene corresponds well with the height of peak solar heating. The RMSE between the pair (1-D, 3-D) of heating rate profiles generally increases with CC and reaches a maximum value of 1.5 K day −1 for overhead sun in the solar spectral range and a maximum value of 3.0 K day −1 in the thermal spectral range (Fig. 7, top middle panel). Although the RMSE is a good measure of an averaged difference between the 1-D and 3-D heating rate profiles, locally the difference between 1-D and 3-D can be much larger. For the stratocumulus scene with CC of 98.9 %, for example, the cloud top cooling is overestimated by about 15 K day −1 and the cloud base warming is overestimated by 5 about 10 K day −1 in the 1-D thermal experiment. The discrepancy between the 1-D and 3-D profile in the thermal spectral range is quantitatively larger than in the solar spectral range and dominates the daytime RMSE at all CC. Nevertheless, during daytime, the difference between the 1-D and 3-D solar heating rate partially compensates the corresponding thermal heating rate difference. The degree of compensation is smallest at SZA of 80 • and increases with increasing solar elevation (Fig. 7, bottom middle panel). The dependence of daytime RMSE on SZA is generally stronger at larger CC.

10
In the solar experiments over the ocean, the 3-D benchmark heating rate within the cloud layer is in general somewhat lower than that in the experiments over land, although this effect is prevailing at small CC and becomes less apparent at larger CC.
At overhead sun, for example, a peak heating rate is reduced by 0.5 K day −1 , 0.6 K day −1 and 0.1 K day −1 in the experiments   In the solar experiments over the ocean (not shown), the 3-D net surface flux is generally larger than that in the experiments over land (e.g., for the cloud scene with CC of 12.0 %, a benchmark value of ∼1040 W m −2 is found at SZA of 0 • and ∼470 W m −2 at SZA of 60 • ). Interestingly, the absolute bias of the 1-D experiment over the ocean at a given SZA and CC is increased (compared to its counterpart over land) by such an amount that the relative bias of the 1-D experiment stays approximately the 10 same as over land (with a deviation of less than 2 %), hinting that the conclusions regarding the surface bias drawn in section 3.3 can be generalized to the entire set of diverse cumulus scenarios.
Finally, the reader should keep in mind, that the results presented in this section should not be interpreted solely as a function of CC. With increasing CC of the scenes from the set, the cloud optical thickness increases as well. This is because both the geometrical thickness and cloud liquid water increase during the evolution of the cloud field (prior to the rain formation), which 15 would also be expected in the real world. Further, apart from CC and LWP, there are plenty of other factors that change from scene to scene and affect the outcome of 3-D experiments and thus the bias of the 1-D calculation. These factors include 3-D cloud geometry, the number of individual clouds in the domain and their spatial distribution.

Statistical synthesis and dependence on cloud liquid water path
In order to obtain a larger data set and thus at least partially overcome the issues discussed in the last paragraph of the previous 20 section, we slightly change the methodology that has been used so far. Namely, for each of the ten cloud scenes, we sample a number of sub-scenes ("windows"), with a horizontal size of 2.8 x 2.8 km 2 at various (x, y)-coordinates within the domain.
In order to create the 1-D experiment, the cloud optical properties within each window are averaged in the same manner as . It is immediately apparent that there is a strong correlation between CC and LWP. As previously suggested, this means that the findings regarding the dependencies on total cloud cover of selected cloud scenarios might be obtained similarly if the cloud scenarios were represented in terms of their optical thickness (but for the sake of brevity we refer to this dependence solely as "CC dependence").
On the whole, the analysis of multiple windows confirms the conclusions that have been drawn in section 3.4 qualitatively (regarding the bias dependence on SZA, A, CC), but extends the range of bias quantitatively. Thus the RMSE between the pair  Although not immediately apparent from the scatter plots, the mean ICA bias is quantitatively lower than the mean 1-D 25 bias, which suggests that the poor representation of clouds in terms of horizontal structure and vertical overlap affects the 1-D surface bias as well.
Examining the dependence on LWP (at fixed values of other parameters) we find that the 3-D solar net surface flux decreases with increasing LWP (at least for overhead sun). The 3-D thermal net surface flux (at a fixed CC) shows little dependence on LWP. Similarly, the dependence of surface biases on LWP is difficult to elucidate. Further investigation is needed to better 30 quantify these effects.

Summary and conclusion
The interaction between radiation and clouds represents a source of uncertainty in numerical weather prediction (NWP) due to both intrinsic problems of 1-D radiation schemes and poor representation of clouds. The underlying question addressed in this study is how large is the bias of radiative heating rates and surface fluxes in NWP models for shallow cumulus clouds and how does it scale with various input parameters of radiation schemes, such as solar zenith angle (SZA), surface albedo (A), total 5 cloud cover (CC) and cloud liquid water path (LWP).
In order to tackle these queries, a set of radiative transfer calculations was carried out for a realistically evolving LES shallow cumulus cloud field, where cloud cover and cloud optical thickness increase with simulation time. For the study we extracted ten time steps with total cloud cover between ∼10 % and ∼100 %. The benchmark experiment was performed on the LES highly-resolved cloud field using a 3-D Monte Carlo radiation model (abbreviated to "3-D" experiment). In order to mimic 10 poor representation of shallow cumulus in NWP models each cloud field was horizontally averaged over the cloudy part of the boxes with dimensions comparable to NWP horizontal grid spacing (several km) and the common δ-Eddington two-stream method with maximum-random overlap assumption for partial cloudiness was applied (abbreviated to "1-D" experiment). An additional experiment was conducted with the same parameter settings as 3-D, except that the Monte Carlo model was run in independent column mode (abbreviated to "ICA" experiment). In other words, the ICA experiment preserves the LES cloud 15 structure and only misses horizontal photon transport, whereas the 1-D experiment misrepresents real cloud structure and lacks horizontal photon transport as well. The comparison between 1-D and 3-D experiments is used to assess the overall bias of NWP radiative quantities (focus of this study), while the comparison between ICA and 3-D experiments allows to separate the effects of horizontal photon transport from those of cloud structure. Each trio (3-D, 1-D, ICA) of experiments was performed in both thermal and solar spectral range. In addition, SZA was varied from 0 • to 80 • with a step of 10 • and different values of 20 shortwave surface albedo (land, ocean) were used.
The vertical profile of the radiative heating rate directly influences atmospheric stratification. Systematic differences in cloud-layer heating rate were found between 1-D and 3-D experiments. In the solar experiments at higher sun elevations (SZA less than 60 • , although this depends slightly on CC) as well as in the thermal experiment, the bias of the 1-D profile shows pronounced vertical gradient within the cloud layer and changes its sign approximately at the height of maximal cloud fraction. 25 In the top part of the cloud layer (i.e., above maximal cloud fraction) the 1-D solar heating rate is too high, while in the bottom part of the cloud layer it is too low, compared to its 3-D counterpart. In the thermal spectral range the opposite is the case, but the effect is quantitatively larger and dominates the total effect of solar and thermal spectral range (at all SZAs). Thus, during nighttime and daytime, the bias of the 1-D heating rate enhances destabilization of the cloud layer by an overestimation of the cooling at cloud top and an overestimation of the warming at cloud bottom (a difference of about −15 K day −1 between 30 1-D and 3-D is observed locally for stratocumulus scenarios). Interestingly, the systematic difference between the 1-D and 3-D solar heating rate is practically insensitive to the choice of surface albedo (i.e., land versus oceanic albedo). In addition to atmospheric heating rates, net surface radiative flux has been investigated with the outcome that the 1-D solar experiment overestimates the corresponding 3-D experiment at low sun (bias up to 80 % over land and ocean), while at high sun the opposite is the case (bias up to −40 % over land and ocean). This positive bias at low sun is largest at intermediate CC, while negative bias at high sun is peaked at larger CC (80 % and beyond). In the thermal spectral range, the 1-D experiment overestimates the amount of 3-D surface cooling, with the maximal bias of about 50 % peaked at intermediate CC.
Overall, the ICA experiment performs better than the 1-D experiment (with respect to the same benchmark). For the abovementioned case of stratocumulus scenarios, for example, a maximum difference of less than 1.5 K day −1 between ICA and 5 3-D is observed locally within the cloud layer. Therefore, one can conclude that resolving horizontally heterogeneous clouds leads to more accurate radiative heating rates than using overlapping fractional plane-parallel clouds in NWP-sized columns.
Since there is a long way to go before shallow cumulus clouds will be resolved within NWP, the aforementioned conclusion implies that the current development of NWP radiation schemes should go hand in hand with the development of advanced cloud schemes generating subgrid-scale cloud structure as realistically as possible. This result is consistent with the findings of 10 the third phase of the Intercomparison of Radiation Codes in Climate Models (ICRCCM III, Barker et al., 2003), in which ICA models outperformed the plane-parallel cloud-overlap models. The work of ICRCCM III assessing only solar radiative transfer was later extended to the thermal part of the spectrum by Kablick III et al. (2011), which brought the same conclusions. Taken together, the results of Barker et al. (2003), Kablick III et al. (2011) and the present study hint that among most promising one-dimensional radiation schemes could be the McICA algorithm (Barker et al., 2002;Pincus et al., 2003), the Tripleclouds 15 method (Shonk and Hogan, 2008) or any other method accounting for unresolved cloud variability. The Tripleclouds method, for example, is an approach to better represent cloud horizontal inhomogeneity by using two regions in each model grid box to represent the cloud as opposed to one. One of these regions represents the optically thicker part of the cloud and the other represents the optically thinner part. The novel optimizations of the Tripleclouds method are currently being developed by the corresponding author of this manuscript and will be discussed in a subsequent study.

20
The question that needs to be addressed next is to what extent do our findings for shallow water clouds apply to other cloud types. A meaningful extension of the present study would include the analogous analysis of ice clouds and mixed-phase clouds.
Especially multi-layered cloud systems which form in the environment with strong vertical wind shear, where the maximumrandom overlap assumption is supposed to break down, appear to be a fruitful avenue for future research. Di Giuseppe and Tompkins (2005) studied the impact of cloud cover on solar radiative biases in deep convective regimes and showed that even 25 apparently complex 3-D convective cloud scenes can often be considered simply in terms of their quasi-two-dimensional cirrus anvil deck, which is an encouraging result.
Nevertheless, full solution for the multiple issues of radiation schemes and their cloud-related problems in today's weather (and climate) models remains a demanding task. This is especially true at the resolution of today's regional (limited-area) numerical models, where a potential 3-D radiation parameterization should take both grid-scale and subgrid-scale radiative 30 effects into account. This is beyond the scope of the present study, but should be perceived as a stimulator for further research on radiation-cloud interactions.
of the manuscript and Dr. Fabian Jakub for providing us with the UCLA-LES cloud field data set. Special thanks to Prof. George Craig for his generous support and fruitful discussion at the initial phase of this work.