The incorporation of the Tripleclouds concept into the δ -Eddington two-stream radiation scheme: solver characterization and its application to shallow cumulus clouds

. The treatment of unresolved cloud-radiation interactions in weather and climate models has considerably improved over the recent years, compared to conventional plane-parallel radiation schemes, which previously persisted in these models for multiple decades. One such improvement is the state-of-the-art Tripleclouds radiative solver, which has two cloudy and one cloud-free region in each vertical model layer and is thereby capable of representing cloud horizontal inhomogeneity. Inspired by the Tripleclouds concept, primarily introduced by Shonk and Hogan (2008), we incorporated a second cloudy region into 5 the widely employed δ -Eddington two-stream method with maximum-random overlap assumption for partial cloudiness. The inclusion of another cloudy region in the two-stream framework required an extension of vertical overlap rules. While retaining the maximum-random overlap for the entire layer cloudiness, we additionally assumed the maximum overlap of optically thicker cloudy regions in pairs of adjacent layers. This extended overlap formulation implicitly places the optically thicker region towards the interior of the cloud, which is in agreement with the core-shell model for convective clouds. The method 10 was initially applied on a shallow cumulus cloud ﬁeld, evaluated against a three-dimensional benchmark radiation computation. Different approaches were used to generate a pair of cloud condensates characterizing the two cloudy regions, testing various condensate distribution assumptions along with global cloud variability estimate. Regardless of the exact condensate setup, the radiative bias in the vast majority of Tripleclouds conﬁgurations was considerably reduced compared to the conventional plane-parallel calculation. Whereas previous studies employing the Tripleclouds concept focused on researching the top-of- 15 the-atmosphere radiation budget, the present work pioneeringly applies the Tripleclouds to atmospheric heating rate and net surface ﬂux. The Tripleclouds scheme was implemented in the comprehensive libRadtran radiative transfer package and can be utilized to further address key scientiﬁc issues related to unresolved cloud-radiation interplay in coarse-resolution atmospheric models.

1 Introduction 20 Radiation schemes in coarse-resolution numerical weather prediction and climate models, commonly referred to as general circulation models (GCMs), have traditionally been claimed to be impaired by the poor representation of clouds (Randall et al., 1984(Randall et al., , 2003(Randall et al., , 2007. Undoubtedly, one of the most rigorous assumptions that persisted in GCMs for multiple decades, was to represent the optically thicker part of the cloud and the other represents the remaining optically thinner part − the method therefore captures cloud horizontal inhomogeneity. Together with the cloud-free region, the radiation scheme thus has three 60 regions at each height and is referred to as the "Tripleclouds" (TC). In the pioneer work of SH08 the layer cloudiness was split into two equally-sized regions and the corresponding pair of cloud condensates (e.g., liquid water content, LWC) was generated on the basis of known LWC distribution. The method was initially tested on high-resolution radar data, where the exact position of the three regions was passed to the radiative solver, capable of representing an arbitrary vertical overlap. In practice, a host GCM usually provides only mean LWC and no information about vertical cloud arrangement. In order to make the method applicable to GCMs, Shonk et al. (2010) derived a global estimate of cloud horizontal variability in terms of fractional standard deviation (FSD), which can be used to split the mean LWC into two components along with the LWC distribution assumption.
Further, they incorporated a generalized vertical overlap parameterization, called the exponential-random overlap, accounting for the aforementioned problematics in strongly sheared conditions. Recently, the method was successfully implemented in the ecRad package (Hogan and Bozzo, 2018), the current radiation scheme of ECMWF Integrated Forecast System (IFS). In 70 contrast to the McICA, which is still operational also at ECMWF, the TC scheme does not produce any radiative noise. As suggested by Hogan and Bozzo (2016) this superiority could become even more valuable in the future if an alternative gas-optics model with fewer spectral intervals than the current RRTM-G (Mlawer et al., 1997) will be developed, since this would increase the level of the McICA noise, but it would not affect the Tripleclouds.
Before the TC solver can be operationally employed, however, it has to be further validated. Whereas all previous studies 75 employing the TC scheme examined primarily the top-of-the-atmosphere (TOA) radiation budget, the present work is aimed at evaluating the atmospheric heating rate and net surface flux. To that end, building upon the Tripleclouds idea, we developed the δ-Eddington two-stream method for two cloudy and one cloud-free region at each height (Fig. 1, bottom right). The prime focus of this paper is to document our exertion of the Tripleclouds concept into the two-stream framework as well as the subsequent implementation of the radiative solver in the comprehensive radiative transfer package libRadtran (Mayer and 80 Kylling, 2005;Emde et al., 2016). Another aim of the present study is to explore the TC potential for shallow cumulus clouds, applying various solver configurations diagnosing atmospheric heating rate and net surface flux. The challenge is to optimally set the condensate pair characterizing the two cloudy regions and geometrically split the layer cloudiness. We test the validity of global FSD estimate in conjunction with three assumptions for sub-grid cloud condensate distribution, which is of practical importance for the application in weather and climate models. 85 The manuscript is organized as follows: In Sect. 2 we introduce the shallow cumulus case study as well as preliminary radiative transfer experiments, demonstrating the importance of representing cloud horizontal heterogeneity. In Sect. 3 we present our version of the TC radiation scheme. Sect. 4 revises existing methodologies for generating cloud condensate pair.
The TC performance is evaluated in Sect. 5. A brief summary and concluding remarks are given in Sect. 6.
2 Cloud data and methodology 90 2.1 Shallow cumulus clouds 2.1.1 Core-shell model for convective clouds We provide a brief note regarding the horizontal distribution of cloud condensate in convective cloud systems. This knowledge will be exploited later when constructing the Tripleclouds radiation scheme. Shallow cumulus clouds are convective clouds, which are often treated with the "core-shell model" Heiblum et al., 2019). In this model the convective 95 "cloud core" associated with updraft motion and increased condensate loading is located in the geometrical centre of the cloud, surrounded by the "cloud shell" associated with downdrafts and condensate evaporation. The core-shell model is supported by multiple observational studies (e.g., Heus et al., 2009;Rodts et al., 2003;Wang et al., 2009) and numerical modelling investigations (e.g., Heus and Jonker, 2008;Jonker et al., 2008;Seigel, 2014) and hence represents the essence of several convection parametrizations. Heiblum et al. (2019) showed that the core-shell model is valid for about 90 % of the typical 100 cloud's lifetime, with the largest discrepancy from the assumed core-shell geometry occurring during the dissipation stage of the cloud. Whereas most of the clouds contain a single core, larger clouds can possess multiple cores. Similarly, clouds in a cloud field have multiple cores, whereby their aggregate effect can be modelled with a core-shell model (Heiblum et al., 2019).

Shallow cumulus cloud field case study
Input for radiative transfer calculations is a shallow cumulus cloud field with a total cloud cover of 54.8 % (visualized in Fig.   105 2), simulated with the University of California, Los Angeles large-eddy simulation (UCLA-LES) model (Stevens et al., 2005;Stevens, 2007). The horizontal domain size is 51.2 x 51.2 km 2 , with the vertical extent of the domain being 3.5 km. A constant horizontal grid spacing of 100 m is applied, whereas the vertical grid spacing is variable ranging from 50 m at the ground to 84 m at domain top. Further details about the UCLA-LES setup can be found in Jakub and Mayer (2017). A 3-D LWC distribution was extracted from a simulation snapshot and the corresponding effective radius (R e ) was parameterized according to Bugliaro are located in the interior of individual clouds, which conforms to the core-shell model (see also Fig. 3). Vertical profiles of averaged LWC, its standard deviation (σ LW C ; simplest measure of cloud horizontal inhomogeneity) and cloud fraction are shown on Fig. 2  The radiative transfer experiments were performed using the libRadtran software (www.libradtran.org), which contains several radiation solvers. The benchmark calculations were performed with the 3-D model MYSTIC, the Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres (Mayer, 2009), which can be run in ICA mode as well. Further, we employed the classic δ-Eddington two-stream method (Joseph et al., 1976) suitable for horizontally homogeneous layers 120 (either fully cloudy or fully clear-sky) and the extension of this method, which allows for partial cloudiness. The latter is the δ-Eddington two-stream method with maximum-random overlap assumption, which was recently implemented into libRadtran in the configuration as described inČrnivec and Mayer (2019) and is ideally siuted as a proxy for the conventional GCM radiative solver (additional instructions regarding its usage are provided in Appendix A).

125
The background thermodynamic state was the US standard atmosphere (Anderson et al., 1986). The parameterization of Hu and Stamnes (1993) was used to convert LWC and R e into cloud optical properties. The solar experiments were performed for a solar zenith angle (SZA) of 0 • , 30 • and 60 • and a surface albedo of 0.25. In the thermal part of the spectrum the surface was assumed to be nonreflective. The shortwave calculations applied 32 spectral bands of the correlated k-distribution by Kato et al. (1999), whereas the longwave calculations employed 12 spectral bands adopted from Fu and Liou (1992). In the Monte Carlo 130 experiments the standard forward and the efficient backward photon tracing were employed in the solar and thermal spectral range respectively. The resulting Monte Carlo noise of domain-averaged quantities is negligible (less than 0.1 %).

Diagnostics and error calculation
The radiative diagnostics include atmospheric heating rate and net (difference between downward and upward) surface flux.
Each diagnostic was examined in the solar, thermal (nighttime effect) and total (daytime effect) spectral range. The error is 135 given by the absolute bias (Eq. 1), relative bias (Eq. 2) and for the atmospheric heating rate additionally by the root mean square error evaluated throughout the vertical extent of the cloud layer (Eq. 3): where y represents the biased quantity and x represents the benchmark.

Preliminary radiative transfer experiments
We present a set of preliminary radiative transfer experiments (listed in Table 1), introducing the 3-D benchmark, the ICA and 145 the conventional GCM calculation. Further, we aim to quantify the various error sources of GCM radiative heating rates, in particular the error related to neglected cloud horizontal heterogeneity.

Benchmark heating rate
The benchmark calculation using MYSTIC (abbreviated to "3-D" experiment) was performed on the highly-resolved LES cloud field (Fig. 4, left). Supposing that the entire LES domain is contained within one GCM column, the quantity of interest 150 is a single vertical profile of radiative heating rate, thus results were horizontally averaged across the domain. Figure 5 (left) shows the resulting benchmark profiles.
In the solar experiment for overhead Sun (Fig. 5, top left) there is a large absorption of radiation in the cloud layer, resulting in a peak heating rate of 10.8 K day −1 . The latter is reached at a height of 1.6 km, which is slightly above the height of maximal cloud fraction (Fig. 2, right). With decreasing Sun elevation the solar heating rate diminishes, exhibiting the maximum of 9.4 155 K day −1 and 5.5 K day −1 at SZA of 30 • and 60 • , respectively. The height where the peak heating is reached stays the same at all SZAs. In the thermal spectral range (Fig. 5, bottom left) the cloud layer is subjected to strong cooling, reaching a peak value of 17.7 K day −1 attained at the same height as the maximum solar heating. Below this height, the magnitude of cooling decreases towards the cloud base, where a slight warming effect is observed.

160
In order to mimic the conditions in conventional GCM models (Fig. 4, right), the cloud optical properties in each vertical layer were horizontally averaged over the cloudy part of the domain, creating a suite of plane-parallel partially cloudy layers.
Consequently, the δ-Eddington two-stream method with maximum-random overlap assumption was employed (abbreviated to "GCM" experiment).
The main shortcomings of the GCM compared to the benchmark (Fig. 5, right) are as follows. In the solar spectral range the 165 peak heating rate is overestimated by 2.7, 2.1 and 0.8 K day −1 at SZA of 0 • , 30 • and 60 • , respectively. In the thermal spectral range the GCM bias artificially enhances radiatively driven destabilization of the cloud layer by an overestimation of cooling by 6.0 K day −1 at cloud-layer top and an overestimation of warming by 3.4 K day −1 at cloud-layer bottom. The GCM error sources are multiple: the misrepresentation of realistic cloud structure, the neglected sub-grid horizontal photon transport as well as the intrinsic difference between the Monte Carlo and two-stream radiative solvers.

ICA and its limitations
To quantify the effect of neglected horizontal photon transport, we run the Monte Carlo radiative model in independent column mode on the original cloud field preserving its LES resolution (Fig. 4,left), with the result horizontally averaged over the domain (abbreviated to "ICA" experiment). Similarly, we applied the δ-Eddington two-stream method within each independent column of the original LES grid (Fig. 4, left) and subsequently averaged the result horizontally (abbreviated to "TSM" 175 experiment). The difference between the ICA and 3-D is a measure of horizontal photon transport. The difference between the TSM and 3-D is a measure of both the horizontal photon transport as well as the intrinsic difference between the Monte Carlo and two-stream radiative solvers.
As anticipated, both indepedent column experiments (ICA, TSM) perform similarly (Fig. 5,right), implying that the intrinsic difference between the radiative solvers is small. Therefore only the ICA is discussed hereafter. The solar bias increases with 180 descending Sun (cloud side illumination; Hogan and Shonk, 2013;Mayer, 2015, 2016), reaching a maximum of −0.7 K day −1 at SZA of 60 • . The amount of thermal cooling is underestimated in the ICA (up to 1 K day −1 ), since realistic cloud side cooling is neglected Mayer, 2014, 2016). Nevertheless, the ICA still overall performs considerably better than the conventional GCM, implying that the major error source of GCM heating rate stemms from the misrepresentation of cloud structure, and not from the neglected horizontal photon transport.

Cloud horizontal heterogeneity effect
In order to isolate the effects of neglected cloud horizontal heterogeneity in a conventional GCM from other effects related to the misrepresentation of cloud structure (e.g., vertical overlap assumption), we employed the GCM radiative solver on the cloud field preserving its LES resolution, but with removed horizontal heterogeneity (Fig. 4, middle). In this way the averaged (plane-parallel) cloud optical properties were applied in each vertical layer, but the realistic 3-D cloud field geometry was 190 retained. The results were horizontally averaged (abbreviated to "HOM" experiment).
The radiative heating rate in the HOM experiment ( can be removed with Tripleclouds? In other words, how well can the continuous probability density function (PDF) of layer LWC be represented by just two cloudy regions (a two-point PDF)?
3 The Tripleclouds radiative solver We first explain in Sect. 3.1 the δ-Eddington two-stream radiation scheme and introduce the terminology. We focus only on those aspects of the method, important to understand its extension to multiple (three) regions, explained in subsequent Sect.
3.2 and 3.3. Differences between the radiative solver of SH08 and our implementation are summarized in Sect. 3.4. Further 200 technical instructions regarding the Tripleclouds usage within the scope of libRadtran are provided in Appendix A.

δ-Eddington two-stream method
In the classic two-stream approach, the entire radiative field is approximated solely with direct solar beam (S) and two streams of diffuse radiation: the downward (E ↓ ) and upward (E ↑ ) component. The widely employed δ-Eddington approximation is a reliable way to account for a strong forward-scattering peak of cloud droplets (Joseph et al., 1976;King and Harshvardhan, 205 1986; Stephens et al., 2001). For the calculations in a vertically inhomogeneous atmosphere, the atmosphere is divided into a number of homogeneous layers, each characterized by its set of constant optical properties. Considering a single layer (j) located between levels (i-1) and (i) (illustrated in Fig. 6) 1 , a system of linear equations determining the fluxes emanating from the layer as a function of fluxes entering the layer can be written as: a 11 a 12 a 13 a 12 a 11 a 23 0 0 a 33 The coefficients a kl in Eq. (4) are referred to as Eddington coefficients. They depend on the optical properties of layer (j) and have the following physical meaning: a 11 -transmission coefficient for diffuse radiation, a 12 -reflection coefficient for diffuse radiation, a 13 -reflection coefficient for the primary scattered solar radiation,

215
a 23 -transmission coefficient for the primary scattered solar radiation, a 33 -transmission coefficient for the direct solar radiation. 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, while index j is used for layer variables. The N vertical layers, enumerated from 1 to N , are enclosed by (N +1) vertical levels, enumerated from 0 to N . For the inclusion of thermal radiation in Eq. (4), the reader is referred to Zdunkowski et al. (2007). The individual layers are coupled vertically by imposing flux continuity at each level. Taking the boundary conditions at TOA (Eq. 5) and at the ground (Eq. 6, with A g representing ground albedo) into account, the radiative fluxes throughout the atmosphere are computed by solving the matrix problem (Coakley and Chylek, 1975;Wiscombe and Grams, 1976;Meador and Weaver, 1980;Ritter and Geleyn, 1992). Henceforth, the calculation of heating rates 225 is straightforward.

δ-Eddington two-stream method for three regions at each height
Consider now a model layer located between levels (i-1) and (i) divided into three regions (Fig. 7). Such layer is characterized by three sets of optical properties and corresponding Eddington coefficients: one for the region of optically thick cloud (superscript "ck"), the other for the region of optically thin cloud (superscript "cn") and the third for the cloud-free region (superscript 230 "f "). In order to apply vertical overlap rules the radiative fluxes corresponding to each of the three regions need to be defined separately at each level (e.g., S ck , S cn and S f ; and analogously for both diffuse components). Total radiative flux at level (i) is thus the sum of both cloudy and the cloud-free component: Figure 7. A model layer between levels (i-1) and (i) divided into three regions.
The equation system (4) is replaced by: 245 so that the fluxes emanating from a certain region of the layer under consideration (e.g., region of optically thick cloud) generally depend on a linear combination of the incoming fluxes stemming from each of the three regions in adjacent layers.
The coefficients starting with T appearing in Eqs. (10), (11), (12) are referred to as the overlap (transfer) coefficients and correspond to layer (j). The coefficient T ck,cn ↓ (j), for example, represents the fraction of downward radiation that leaves the base of optically thick cloud of layer (j-1) and enters the optically thin cloud of layer under consideration (j). The overlap 250 coefficients quantitatively depend on the choice of the overlap rule, which will be discussed in the next subsection (3.3). For a three-region layer, the boundary condition at TOA (Eq. 5) implies: The boundary condition at the ground (Eq. 6) is extended to: which assumes that the downward fluxes leaving the lowest model layer, after reflection enter the same sections of individual cloudy and cloud-free air (isotropic ground reflection).

Overlap considerations
The layer cloud fraction C is given by: In our implementation we demand the following relationship between the individual cloud fraction components: where α is a constant between 0 and 1. We apply the widely used maximum-random overlap assumption (Geleyn and Hollingsworth, 1979) for the entire layer cloudiness (sum of optically thick and thin cloudy regions), where adjacent cloudy layers exhibit maximal overlap and cloudy layers separated by at least one cloud-free layer exhibit random overlap. If the 275 cloudy layers are splitted into two parts, however, this overlap rule is not sufficient and needs to be extended. Therefore, we additionally assume the maximum overlap of adjacent optically thicker cloudy regions and abbreviate this extended overlap rule to the "maximum 2 -random overlap". This assumption implicitly places the optically thicker cloudy region towards the interior of the cloud in the horizontal plane, which is in line with the core-shell model. Now we can quantitatively determine the overlap coefficients in Eqs. (10), (11) and (12) for the maximum 2 -random overlap. 280 We consider the transmission of downward radiation through two adjacent layers with partial cloudiness. Four possible geometries, illustrated in Fig. 8, need to be treated. For the situation depicted on the top left panel of Fig. 8, the transmission of direct radiation can be formulated as follows. The optically thick cloud of layer (j-1) transmits S ck (i-1), the optically thin cloud transmits S cn (i-1) and the cloud-free region transmits S f (i-1). These three components of the transmitted radiation must then 13 https://doi.org/10.5194/acp-2020-334 Preprint. Discussion started: 30 April 2020 c Author(s) 2020. CC BY 4.0 License. be distributed between the three regions of the lower layer (j). The maximum overlap of optically thick cloudy regions implies 285 that the entire radiation S ck leaving the base of layer (j-1) enters the optically thick cloud below: and none of it enters the other two regions:
To ensure the maximum overlap of cloudy layers as a whole, the remaining cloudy flux at the base of layer (j-1), namely the S cn (i-1), needs to be lead into the two cloudy regions of the lower layer, with the priority to enter the optically thick cloud.

295
T cn,cn T cn,f The cloud-free flux S f at the base of layer (j-1) is distributed according to: 305 We leave the derivation of overlap coefficients for other three geometries as an exercise for the reader. The transmission of upward radiation is managed via overlap coefficients T a,b ↑ in a similar fashion. It should be noted that the same coefficients govern the reflection, whereby the upward reflection of downward radiation is treated with T a,b ↓ and the reverse situation is treated with T a,b ↑ . Pairwise overlap as employed here ensures that the matrix problem is fast to solve. Whereas a drawback of the core-shell model and thereby the outlined overlap is that it underperforms in case of vertically developed cloud systems in 310 strongly sheared conditions, the present Tripleclouds implementation is an excellent tool to study shallow convective clouds.
In this way the effects of cloud horizontal inhomogeneity are tackled in isolation, while the issues related to vertical shear are eliminated.
The Tripleclouds radiative solver has been successfully implemented in the libRadtran package. Technically, the calculation of overlap coefficients is performed in an autonomous function enabling flexible modifications of overlap rules in the future.

Differences to the radiative solver of SH08
We briefly outline the main differences between our radiative solver and that of SH08 regarding the incorporation of threeregion layers in two-stream equations. SH08 used another version of a two-stream solver implemented in the radiation scheme devised by Edwards and Slingo (1996). The way, how the two-stream solver was incorporated in the Edwards-Slingo code, resulted in a complicated expression for upward fluxes, when the solver was extended to multiple regions (their Eq. 15). There-320 fore SH08 simplified this expression, achieving higher computational efficiency, but bringing certain physical shortcomings. They explained: "Physically, this means that, at each height, the downwelling radiation in a particular region is either reflected back or absorbed; none is reflected up into another region." The Tripleclouds radiative solver as constructed in the present work allows for a full interaction of the three regions, meaning that downward radiation in a particular region being reflected upward (and vice versa), can in general be distributed between 325 all three regions. This feature is illustrated in Fig. 9 and could play a comparatively important role especially in the solar calculations, because of significant scattering effects. Moreover, this physical mechanism was recently recognized as relevant and incorporated also into the latest version of the Tripleclouds solver at ECMWF (Hogan et al., 2019). Whereas the work of Hogan et al. (2019) represents the first study utilizing the complete TC form researching the TOA cloud radiative forcing, the present study pioneeringly applies the full TC scheme to atmospheric heating rate and net surface flux.

Methodologies to generate the LWC pair
In order to apply the TC radiative solver, a pair of LWC characterizing optically thin and thick cloudy regions (LWC cn , LWC ck ) needs to be created in each vertical layer. In Sect. 4.1 we revise the original Tripleclouds method introduced by SH08, later referred to as the "lower percentile method" (Shonk et al., 2010), which can only be applied if the LWC distribution is known.

The lower percentile method
In this method it is assumed that the LWC distribution in each vertical layer can be approximated with the normal distribution: where LW C is layer mean LWC and σ LW C is its standard deviation. The distribution of LWC is divided into two regions through a given percentile of the distribution, denoted as "split percentile (SP)". The latter is chosen to be the 50th percentile or 340 the median, which splits the cloud volume into two equal parts (i.e., cloud fraction in each vertical layer is halved). The LWC of the optically thin cloud (LWC cn ) is determined as the value corresponding to the so-called "lower percentile (LP)" of the distribution. This is chosen to be the 16th percentile based on the following considerations. We adjust the two LWC values in a way that the mean LWC in the layer is conserved: 345 and that they are separated by two standard deviations: For a Gaussian distribution, the latter constraint has a desired property that the variability within each of the two cloudy regions (measured by σ LW C ) is the same as that within the entire cloud in the layer. Equations (32) and (33) give the following relationship for LWC cn : The fraction of the distribution with LWC lower than LWC cn is therefore: which corresponds to the LP of 16. Finally, the LWC ck is determined using Eq. (32) to conserve the mean. Figure 10 shows the resulting LWC pair when the LP method is applied on shallow cumulus cloud field.

355
It should be noted that the choice of the 16th percentile as the LP and the 50th percentile as the SP is based solely on theoretical considerations. In practice, the LP and SP are the two tunable parameters, that can be adjusted according to their performance on real cloud data. Even though the optimal setting varies, SH08 exposed that the combination of LP of 16 and SP of 50 generally serves well in both solar and thermal spectral range for vast ranges of cloud data.

360
This method in its initial formulation by Shonk et al. (2010) implicitly assumes that LWC is normally distributed as well.
Thereby the cloudiness in each vertical layer is partitioned into two regions of equal size and the pair of LWC (LWC cn , Figure 10. LWC profiles obtained with the LP method. Figure 11. The actual FSD of the shallow cumulus. The grey-shaded area represents the uncertainty of global FSD estimate, centered around its mean value (black line).
LWC ck ) is obtained by: where FSD represents the fractional standard deviation of LWC: Since in practice only LW C is known within a GCM grid box, the FSD has to be parameterized. A review of numerous studies (Cahalan et al., 1994a;Barker et al., 1996;Pincus et al., 1999;Smith and DelGenio, 2001;Rossow et al., 2002;Hogan and Illingworth, 2003;Oreopoulos and Cahalan, 2005; SH08) carried out by Shonk et al. (2010) gave a globally representative FSD of 0.75 ± 0.18. Figure 11 shows the actual FSD for the present shallow cumulus: although this FSD is strongly dependent 370 on the position within the cloud layer, it predominantly lies within the range of global estimate.
If the cloud condensate is normally distributed, substracting σ LW C from the LW C to obtain the LWC cn in Eq. (36) corresponds approximately with the 16th percentile. For more realistic lognormal and gamma distributions, the 16th percentile (advocated by SH08) is given by relationships presented in Hogan et al. ( , 2019, whereby the LWC ck is again obtained by conserving the layer mean.

375
In order to test the validty of global FSD estimate, we applied its mean value (0.75) to create the pair of LWC in each vertical layer containing cloud. Further, to test the sensitivity of TC radiative quantities to the assumed form of the sub-grid cloud condensate distribution, we employed the FSD method in conjunction with all three distributions (the resulting profiles are shown in Fig. 12).

380
We evaluated the TC radiative solver with both LP and FSD methods. The effective radii characterizing the two cloudy regions were kept the same (averaged R e ). The setup of radiation calculations was as described in Sect. 2.2. The results of the various TC experiments are compared with the conventional GCM, which approximates the cloud condensate distribution with a onepoint PDF and can be perceived as an upper bound for the tolerable TC error. In addition, the ICA, which resolves the full sub-grid PDF, is shown as well. In Sect. 5.1 we examine atmospheric heating rate, whereas in Sect. 5.2 we discuss net surface 385 flux.
5.1 Atmospheric heating rate

Tripleclouds with LP method
We evaluate first the TC radiative solver when the LP method is used to obtain the pair of LWC. The results of this experiment, denoted as "TC(LP)", are shown in Figs. 13 and 14. It is apparent that the TC(LP) is overall significantly more accurate than 390 the GCM. In the solar spectral range for overhead Sun (Fig. 13, top middle), the maximal bias within the cloud layer is reduced from 2.7 K day −1 to only 0.7 K day −1 . Whereas the largest bias reduction is observed within the cloud layer, the heating rate above and below the cloud layer is considerably improved as well, explained as follows. The non-homogeneous clouds have lower mean shortwave albedo and absorptivity than the corresponding plane-parallel cloudiness with the same mean optical depth ( Fig. 2 of Cairns et al., 2000). This implies that the non-homogeneous cloud in the TC configuration reflects less of the 395 incoming solar radiation upward (leading to a reduction of the positive GCM bias above the cloud layer) and simultaneously absorbs less radiation (leading to a reduction of the positive GCM bias in the cloud layer), compared to the homogeneous cloud in the GCM. Consequently, more radiation is transmitted through the cloud layer and absorbed in the region below the cloud layer in the TC experiment compared to that in the GCM, which reduces the negative GCM bias in this region. At SZA of 30 • the behaviour is qualitatively similar, with the maximal bias of 2.1 K day −1 within the cloud layer reduced by a factor of 5.

400
At SZA of 60 • , the maximal bias of 0.8 K day −1 within the cloud layer becomes of the opposite sign, but is still smaller in magnitude (−0.4 K day −1 ), when the TC(LP) is applied in place of the conventional GCM. In the layer above and especially below the cloud layer, however, the bias is slightly increased. Noteworthy, at all three SZAs, the 3-D radiation feature at cloud base (increased heating due to surface reflection of radiation) can not be properly accounted for using the TC solver.
In the thermal spectral range (Fig. 13, bottom middle), the degree of artificially enhanced destabilization of the cloud layer, 405 arising from the overestimation of cloud top cooling and cloud base warming in the GCM, is drastically reduced when the TC(LP) is applied, interpreted as follows. The non-homogeneous clouds have lower mean longwave emissivity and absorptivity than the corresponding homogeneous clouds with the same mean optical depth. Thus the non-homogeneous cloud top in the TC experiment emits less radiation compared to the homogeneous cloud top in the GCM configuration, which reduces the negative GCM bias at cloud top. Similarly, the non-homogeneous cloud base in the TC experiment absorbs less of the radiation stemming 410 from the warmer atmospheric layers underneath the cloud, compared to the homogeneous cloud base in the conventional GCM, which reduces the positive GCM bias at cloud base. As anticipated, in the region above and below the cloud layer, the difference between the TC and the GCM is only marginal.

Tripleclouds with FSD method
We first examine the TC(FSD) experiment when the gaussianity of cloud condensate is assumed − this experiment is consid-415 erably more accurate than the conventional GCM as well. As an illustration, the daytime cloud-layer RMSE of 1.7 K day −1 is reduced to 0.3 K day −1 at SZA of 60 • (Fig. 14). Furthermore, the TC(FSD) experiment is even slightly more accurate than the TC(LP) especially in the thermal spectral range and in the solar spectral range at SZA of 30 • and 60 • , whereas at SZA of 0 • the situation is reversed (Fig. 13, middle column). The largest discrepancy between the two TC experiments is observed in the  When examining the entire set of TC(FSD) experiments it is apparent that the radiative heating rate is considerably more accurate compared to the conventional GCM regardless of the exact assumption for the LWC distribution. Although the Gaussian distribution was ranked worst when fitted to the actual PDF, the gaussianity assumption with global FSD performed best 430 in practice, contemplated as follows. In the central part of the cloud layer around maximum cloud fraction the actual FSD of the present shallow cumulus (0.95) is larger than the assumed global estimate. The latter is primarily due to great amount of cloud side area in this region, an essential characteristic of broken cloud field, which generally contributes to increased variability (Hill et al., 2012(Hill et al., , 2015. Since the assumption of gaussianity implies the largest difference between the LWC pair characterizing the two cloudy regions (Fig. 12), it partially accounts for the missing variability provided by the global estimate.

435
More sophisticated FSD parameterizations are tempting to be tested in the future. bias. This is attributed to the fact that the plane-parallel GCM cloudiness leads to an increased solar absorption and hence reduced cloud-layer transmittance. The latter reduces downward flux reaching the surface and profoundly underestimates the net flux. At nighttime, the plane-parallel cloud in the GCM emits a greater amount of radiation towards the surface compared to heterogeneous cloud in the ICA, leading to a reduction of surface net flux bias.

Net surface flux
When the Tripleclouds is applied either with the LP or the FSD method instead of conventional GCM radiation scheme, the 455 daytime net surface flux bias of −55 W m −2 (or −8 %) is substantially reduced to −5 W m −2 (or −1 %) at overhead Sun and similarly for SZA of 30 • (assuming gaussianity of cloud condensate). At SZA of 60 • and especially at nighttime, radiative bias in the various TC experiments increases compared to the GCM bias. This indicates that the TC in its current configuration should be taken with caution when applied to surface thermal flux, as its usage can lead to degradation of the nocturnal surface budget compared to simple plane-parallel model.

Summary and conclusions
Inspired by the Tripleclouds concept of Shonk and Hogan (2008), we incorporated a second cloudy region in the widely used δ-Eddington two-stream method with maximum-random overlap assumption for partial cloudiness. The resulting radiation scheme thus has two cloudy and one cloud-free region in each vertical layer and is capable of representing cloud horizontal variability. The inclusion of a second cloudy region into the two-stream framework required an extension of vertical overlap 465 rules. While retaining the maximum-random overlap for the entire layer cloudiness, we additionally assumed the maximum overlap of optically thicker cloudy regions in pairs of adjacent layers. This implicitly places the optically thicker region towards the interior of the cloud in the horizontal plane, while the optically thinner region resides at cloud periperhy, which is in line with the core-shell model for convective clouds.
The constructed Tripleclouds radiative solver was evaluated on a shallow cumulus cloud field. The validity of global esti-470 mate of fractional standard deviation − a common measure of cloud horizontal variability − was tested along with different assumptions for sub-grid cloud condensate distribution (Gaussian, gamma, lognormal), which are frequently applied when parameterizing clouds in weather and climate models. In the vast majority of experiments the Tripleclouds performed better than the conventional plane-parallel GCM scheme. The error of atmospheric heating rate was substantially reduced at daytime and nighttime (up to fivefold cloud-layer RMSE reduction). In case of net surface flux the daytime bias was generally 475 depleted as well, whereas the nighttime bias reduction was less pronounced, suggesting that the computationally more efficient plane-parallel scheme could be retained in this case.
The question that needs to be addressed next is to what extent do our findings for a shallow cumulus case study with intermediate cloud cover apply to a larger set of scenarios comprising a wide range of cloud cover. This question is relevant, because horizontal variability might essentially depend on cloud fraction. Similarly, the degree of cloud horizontal variability 480 might depend on the GCM grid resolution, which has to be investigated in more detail in the future. Furthermore, organizational aspects of shallow convection should be addressed in the context of the present study. Mesoscale shallow convection sometimes occurs in the form of uniformly scattered cumuli, but is also frequently organized into cloud streets, clusters or mesoscale arcs (Agee et al., 1973;Atkinson and Zhang, 1996;Wood and Hartmann, 2006;Seifert and Heus, 2013). The robustness of the results on the nature of cloud organization should be examined next. Recently, Stevens et al. (2019) proposed four mesoscale 485 cloud patterns frequently observed in trade wind regions, which they labeled Sugar, Flower, Fish and Gravel. A follow-up study of Rasp et al. (2019) proved that the four patterns correspond to physically meaningful cloud regimes, each of them being associated with specific large-scale environmental conditions. These climatologically distinct environments should exhibit a highly variable cloud water variance. If this proves true and if the internal cloud variability is properly quantified, a regime-several radiative transfer equation (RTE) solvers. To promote the usage of both recently implemented two-stream solvers (termed "twomaxrnd" and "twomaxrnd3C"), which are both coded in C programming language, basic guidelines are given below. For a complete description on how to set up the background atmosphere and other input parameters, the reader is referred to the libRadtran user manual, which is included in the software package. The output quantities of both algorithms include either radiative fluxes (default) [W m −2 ] or heating rates [K day −1 ]. Whereas examples provided below illustrate the 505 treatment of water clouds, both RTE solvers can be applied to ice clouds in a similar fashion.

A1 RTE solver: "twomaxrnd"
The δ-Eddington two-stream method with maximum-random overlap assumption for partial cloudiness in the configuration as documented in Sect. 2.2 ofČrnivec and Mayer (2019)  where cf.dat is the standard libRadtran file containing cloud fraction vertical profile and wc.dat is the standard 1-D file 515 defining water cloud properties.
A2 RTE solver: "twomaxrnd3C" The Tripleclouds radiative solver, effectively the δ-Eddington two-stream method for two cloudy and one cloud-free region at each height with maximum 2 -random overlap assumption, as described in Sect. 3 of the present work, is invoked as follows: where cf.dat is again the standard file containing the vertical profile of cloud fraction. It is important to note that this file determines the cloud fraction of the entire layer cloudiness (sum of optically thick and thin cloudy regions). The division of the latter into two components is managed via newly introduced parameter twomaxrnd3C_scale_cf, which corresponds to the parameter α in Eqs. 20 and 21. The split of averaged cloud water properties into two components is not yet automated, rather 530 the user is asked to preprocess both cloud files depending on his/her specific needs. The resulting wck.dat and wcn.dat are 1-D water cloud files, defining properties of optically thick and thin cloudy regions, respectively (note that the option profile_file is solely the generalization of the wc_file command).

Appendix B: Analytical probability density functions
We outline the relationship between LW C, σ LW C , FSD (denoted as f LW C in the following) and the parameters used to 535 describe lognormal and gamma distributions, which were applied to fit the modeled LWC distributions.
A lognormal distribution of LWC is defined as: The parameters of the lognormal distribution, LW C 0 and σ 0 , can be defined in terms of LW C and f LW C in the following fashion: σ 2 0 = ln(f LW C + 1).
A gamma distribution of LWC is defined as: 545