Computation of domain‐average radiative flux profiles using Gaussian quadrature

A method for calculating domain‐average radiative flux profiles, called Gaussian Quadrature Independent Column Approximation (GQ‐ICA), is introduced and assessed using cloud properties retrieved from A‐Train satellite data. This method could be suitable for use in large‐scale atmospheric models. Like the Monte Carlo ICA (McICA), GQ‐ICA uses N stochastically generated subgrid‐scale cloudy columns. The independent variable is the sorted, from smallest to largest, sequence of N sub‐column values of liquid and ice cloud water paths. The integrand is essentially the radiative transfer equation. Accurate GQ integration requires integrands to be relatively smooth functions. Unlike McICA, GQ‐ICA performs full solar and infrared spectral integrations on nG <  < N sub‐columns which are identified by rules governing nG‐node GQ. The nG flux profiles are appropriately weighted and summed to give domain averages. Several sorting procedures were considered, and all results are based on the CCCma radiation algorithm. For solar radiation, 1‐node GQ‐ICA can produce significant bias errors, but its random errors are generally less than McICA's. These biases, however, are almost eliminated by 2‐node GQ‐ICA. For GQ‐ICA to better McICA's random errors for infrared fluxes, at least the 2‐node version is needed. Ultimately, 2‐node GQ‐ICA random errors for net fluxes at surface and top‐of‐atmosphere are typically 30–50% of McICA's. This is partly because solar and infrared solvers operate on the same sub‐columns. GQ‐ICA random errors for atmospheric heating rates are comparable to McICA's even for 3‐node GQ‐ICA. Computational times required for the 2‐ and 3‐node GQ‐ICA are, respectively, ∼180 and ∼230% of McICA's.


INTRODUCTION
Clouds play a large part in maintaining and adjusting to alterations in Earth's radiation budget. As such, it is essential that the radiative properties of cloudy atmospheres be represented well in General Circulation Models (GCMs). Despite recent advances, descriptions of clouds, and related processes, in GCMs are still lacking in some respects. This is to be expected given the computational constraints faced by GCMs which usually require their grid cells to be larger than ∼100 × 100 km. That leaves myriad cloud and radiative processes and fluctuations in need of far-reaching parametrizations; in particular, subgrid-scale interactions between radiative transfer and cloud vertical overlap and horizontal inhomogeneity (e.g. Stephens, 1988;Davis et al., 1990;Barker, 1996;2008;Cairns et al., 2000;Hogan and Illingworth, 2000;Bergman and Rasch, 2002;Mace and Benson-Troth, 2002;Pincus et al., 2005;Tompkins and Di Giuseppe, 2007;Naud et al., 2008).
The subgrid-scale radiation-cloud problem received a decided boost from the Monte Carlo Independent Column Approximation (McICA) method (Barker et al., 2002;Pincus et al., 2003). McICA separates descriptions of the unresolved atmosphere-surface system from the radiative transfer solver, as in a conventional Monte Carlo solver (see several chapters in Marshak and Davis, 2005), making for unbiased estimates of domain-average broadband fluxes, relative to the full ICA, as well as relatively simple and flexible coding. Many modelling centres adopted, and continue to use, McICA (Pincus et al., 2006;Räisänen et al., 2007;2008; Barker et al., 2008;Morcrette et al., 2008;Hill et al., 2011) in conjunction with a stochastic cloud generator (Barker et al., 1999;) that supplies sub-columns based on information from the GCM and assumptions regarding the structure of unresolved clouds (Hogan and Illingworth, 2000;Barker, 2008;Oreopoulos et al., 2012;Zhang et al., 2013).
The manner in which McICA samples physical space as it steps through spectral integration, usually with a correlated k-distribution (CKD) routine, leads to radiative noise that gets conveyed back to the GCM. While several experiments using fixed boundary conditions have shown GCMs to be largely immune to McICA's noise (e.g. Barker et al., 2008), the fact that the short-wave (SW) and long-wave (LW) radiative transfer models operate on distinctly different sampled atmospheres creates unnecessary noise for net fluxes and is plainly unphysical.
This study presents an extension of the regular McICA method in an attempt to: (a) bring consistency to the SW and LW transfer solvers, and (b) reduce McICA's radiative noise. The latter point is potentially important for coupled atmosphere-ocean simulations (cf. Räisänen et al., 2008). The basis of the proposed method rests on Gaussian quadrature (GQ) integration with the independent variable being sorted cloud water path, as provided by the stochastic sub-column generator, and the integrands being the radiative transfer equation. This methodology resembles Shonk and Hogan's (2010) tripleclouds scheme.
Section 2 reviews briefly the McICA method and introduces the GQ-based extension of McICA. Section 3 describes the cloudy atmospheric data, inferred from A-Train satellite measurements, as well as the SW and LW radiative transfer models used for this study. There are various ways one can go about sorting cloud water paths for GQ integration and some are discussed and assessed in both section 4 and an appendix. Results are presented in section 5, and the final section offers a brief summary and conclusion.

ESTIMATING DOMAIN-AVERAGE RADIATIVE FLUXES FOR CLOUDY ATMOSPHERES
The primary purpose of this section is to introduce the use of Gaussian quadrature (GQ) as a means of approximating mean radiative flux profiles for large areas that contain cloud. As the new method is an extension of McICA, McICA is reviewed in the following subsection with its GQ counterpart, denoted hereinafter as GQ-ICA, presented in the second. Like regular McICA, GQ-ICA can be applied to domains of well-resolved clouds, such as inferred from satellite observations or produced by cloud-resolving models. Here, however, application was limited to domains with unresolved clouds, in need of stochastic generation, such as in weather and climate prediction models.

Monte Carlo Independent Column Approximation (McICA)
If a domain consists of explicit sub-columns, such as provided by a cloud-resolving model, or virtual sub-columns that are generated stochastically, such as in a GCM, and the total cloud fraction is C tot , full ICA benchmark, domain-average, broadband fluxes are defined as where N 0 and N c are numbers of clear-and cloudy-sky sub-columns, K is the number of spectral intervals, and F 0 i,k and F c i,k are clear-and cloudy-sky fluxes for the ith sub-column. In this study, correlated k-distribution (CKD) radiation algorithms were used with K points in cumulative probability space (CPS: Li and Barker, 2005). Assuming that all N 0 clear-sky sub-columns are the same, domain-average spectrally integrated clear-sky fluxes are F 0 , with F c being the corresponding cloudy-sky value. This implies that only a single clear-sky radiative calculation is needed. One could save computation time by randomly sampling N ′ c < N c cloudy sub-columns and arriving at an unbiased estimate of whilst incurring a random sampling error. Going further and stipulating that each spectral interval receive a randomly sampled cloudy sub-column yields the McICA method (Barker et al., 2002;Pincus et al., 2003) in which where F c i ′ (k),k indicates that for each k there is a random sub-column of index i ′ (k). If one is dealing with a virtual domain of sub-columns, such as in a GCM, Equation 3 becomes where C ′ tot is computed from the generated sub-columns. In this case, cloudy sub-columns are produced by a stochastic generator (e.g.  in which the true, but operationally impractical, benchmark is still Equation 1 but with N c → ∞ random columns. In operational settings, unbiased random errors for F ′ can be substantial (cf. Barker et al., 2008) given that they arise from both sub-column stochastic sampling and arbitrary association of sub-columns with k. While the former source of noise is obvious, the latter can be appreciated by repeated application of Equation 4 to a sample of K sub-columns, but shuffling the sub-columns each time. While C ′ tot would remain set, the sum in Equation 4 would differ with each shuffle. By extension, spectral integrations over SW and LW spectra experience unique sample atmospheres; even if both integrations were to use the same sample of sub-columns. This unphysical aspect of McICA, including Equations 2 and 3, cannot be circumvented, by definition, and enhances random errors for net fluxes, though techniques have been devised to reduce, but not eliminate, them Hill et al., 2011). 1

Gaussian Quadrature Independent Column Approximation (GQ-ICA)
The objective of the method presented in this subsection is to rectify McICA's unphysical sampling by forcing the SW and LW spectral integrations to experience a common sampled atmosphere. It makes direct use of Gaussian quadrature integration (Abramowitz and Stegun, 1972;Li, 2000) and is referred to, hereinafter, as the Gaussian Quadrature Independent Column Approximation (GQ-ICA).
GQ is a means of evaluating the definite integral of a functionf (x) on x ∈ [a, b] via a small number of evaluations of f (x). While there are many forms of GQ, this study used only Gauss-Legendre quadrature for [0, 1], for which the n G -node approximation is where w l and x l ∈ [0, 1] are well-defined weights and values. This expression is an equality when f (x) is a polynomial of degree ≤2n G − 1, but is generally accurate when f (x) can be approximated well by such a polynomial. The GQ method is unsuitable for functions with singularities or large discontinuities. Begin by producing N sub-columns with Räisänen et al.'s (2004) stochastic cloud generator. As with regular McICA, where N c is the number of sub-columns containing cloud. For each cloudy sub-column compute cloud water path W(i) ; i = 1, … , N c , which is the most important variable, next to C ′ tot , governing radiative fluxes. Now, using 5, domain-average cloudy-sky flux is then where the right-hand side of the top line is simply the ICA evaluated for the sorted sub-columns, the second line is the transformed integral equivalent of the ICA, and the third line indicates summation of n G full-spectrum fluxes using the GQ-specified sub-columns; where it is recognized that dx = i/N c − (i − 1)/N c = 1/N c . As such, both the SW and LW routines operate on the same sub-columns. The obvious potential penalty associated with this method is that a portion of the radiative transfer algorithms get evaluated n G times. In terms of CPU usage, n G should be kept as small as possible while maintaining adequate sampling of unresolved variable cloud. How best to sort W is discussed in section 4 and an appendix. First, however, the following section presents the data and radiative transfer models that were used here to assess GQ-ICA and McICA.

EXPERIMENTAL SET-UP
To verify the GQ-ICA proposed in section 2.2 and section 4, descriptions of cloudy-sky domains were used that derived from cloud properties retrieved from measurements made by CALIPSO's dual-wavelength lidar (Winker et al., 2003) and CloudSat's 94 GHz cloud-profiling radar (Stephens et al., 2002) during August 2007. Layer cloud fractions were derived from CloudSat's 2B-GEOPROF-LIDAR product (Mace et al., 2007), while cloud liquid and ice water contents came from 2B-CWC-RO (Austin 2007). Each layer is 0.24 km thick with the lowest two or three layers being unworkable due to ground-clutter contamination. Columns are ∼1.2 km wide, though they are essentially infinitely wide when used by the one-dimensional (1D) radiative transfer models. CloudSat's ECMWF-AUX files report profiles of based on their LWP (green lines). Dots represent cloud-bearing cells. Reflected SW fluxes at TOA for each sub-column are shown by red lines ( 0 was set to 0.9). GQ integration columns, and their respective fluxes, for n G = 1, 2 and 3 are shown by blue lines. Estimated domain-average reflected SW flux is listed above each plot. The ICA benchmark is 341.6 W/m 2 pressure, temperature, and mixing ratios of water vapour and ozone for each column. See Barker (2008) for additional details.
Assuming that GCM grid-columns are nominally ∼250 km across, CloudSat-CALIPSO data were partitioned into 256 km segments. Only segments with C tot > 0.05 were used. This resulted in 3,026 domains distributed across all latitudes. For each domain, 10 values of cosine of solar zenith angle 0 were determined by selecting random times between sunrise and sunset as a function of latitude and date (Barker et al., 2015). Profiles of cloud fraction, mean and variance of cloud water content were computed for each domain and provided to Räisänen et al.'s (2004) stochastic cloud generator which then produced N = 150 sub-columns per domain per radiation calculation. The generator also used vertically constant effective decorrelation length for cloud fraction L * cf (Barker, 2008) defined as the solution to where C ′ tot is total cloud fraction from the generator given layer cloud fraction profile from CloudSat-CALIPSO and N samples, with C tot being computed directly from the satellite data. Decorrelation lengths for cloud water content were set to L * cf ∕2. The broadband radiative transfer models used throughout this study were from the CCCma-GCM (von Salzen et al., 2013). Gaseous absorption for both the SW and LW is handled by the CKD method (Li and Barker, 2005) using 35 and 46 points in CPS for SW and LW bands, respectively. SW transfer is accounted for by the delta-Eddington approximation (Wiscombe, 1977) while LW follows Li (2002) and . Optical properties for liquid cloud are described by Dobbie et al. (1999) for SW and Lindner and Li (2000) for LW radiation. Ice cloud SW properties are described by Fu (1996) and LW by Fu et al. (1998). For simplicity, effective radius for liquid droplets was set to 10 m, while ice crystal effective geometric length was set to 46 m. All surfaces were assumed to be ocean with broadband albedo and emissivity equal to 0.06 and 1.0, respectively.
All simulations, whose results are presented in the following two sections, employed the full Independent Column Approximation (ICA) method to represent benchmark domain-average flux profiles. They were obtained by applying full-spectral integrations of the radiation transfer models just described to each of the 150 sub-columns per domain and averaging. The experiments were either the McICA or GQ-ICA which operated on the same 150 sub-columns.

APPLICATION OF GQ-ICA
The immediate detail facing GQ-ICA is sorting of sub-columns. To illustrate this, consider a domain with liquid clouds only (ice clouds simply zeroed-out). Figure 1 shows the GQ-ICA procedure with N = 150 stochastic sub-columns based on CloudSat-CALIPSO profiles of cloud fraction, mean and variance of cloud water content, and L * cf . Dots indicate cells that contained cloud. Cloudy columns were sorted such that their liquid water paths (LWP) increase monotonically from left to right (green lines). Corresponding upwelling broadband SW fluxes at top-of-atmosphere (TOA) for each sub-column are represented by red lines. They also increase monotonically for this example. Blue lines indicate GQ nodal positions and fluxes for n G = 1, 2 and 3. As expected, and in general, as n G increases, differences between estimates and ICA benchmarks of F c diminish and almost vanish by n G = 3; but they do not necessarily approach zero because the flux curve only approximates a smooth polynomial.
When ice and liquid clouds are present in a domain, the sorting process can become complicated. This is because LWP is usually much greater than ice water path (IWP). Consider now the same domain as used to produce Figure 1, but this time include the upper-level ice cloud. Again, N = 150 stochastic sub-columns were computed resulting in C ′ tot ≈ 0.92. Figure 2a shows the cells that contained cloud as they came out of the generator along with the 12 cloudless sub-columns on the left.
The simplest way to deal with a combination of liquid and ice clouds is to sort LWP + IWP. This approach, referred to as combined sorting, and others are discussed in the appendix and shown to have some limitations. The main issue with combined sorting is that often columns with only ice occur . The upshot is that with so few GQ samples, many ice clouds simply get missed for they lie to the left of the first GQ point. To make this approach even moderately accurate, at least two applications of a 1-node GQ are required; one for ice-only columns and another for the remainder.
Commonly, including the example used for Figure 2, boundary-layer and cirrus clouds have little physical connection (e.g. Cotton and Anthes, 2013). Therefore, the portion of ice cloud overlapping liquid cloud could possibly be sorted on IWP while independently sorting LWP. Figure 2b shows the result of this sorting process which is referred to as double sorting. This produces two regions of sorted ice with the thinner region being operated on by 1-node GQ and the remainder by n G -node GQ. Of course, doing this breaks the vertical correlation as dictated to, and produced by, the stochastic generator. If the portion of ice cloud is very small, however, the randomly realized distribution of IWP is used. The criterion used here to invoke double sorting was that the fraction of a domain with overlapping ice cloud had to exceed 80% of the liquid cloud. If the ice cloud fraction is too small, one of the GQ points in 2-or 3-node GQ integration will miss sub-columns of ice. In terms of radiative performance, the appendix shows that of the four sorting methods tried, SW and LW fluxes are modelled best by double sorting. As such, only results for double sorting are shown in the following section.

RESULTS
In this section the performances of 1-, 2-and 3-node renditions of GQ-ICA, using only double sorting, were assessed along with the regular McICA for the dataset described in section 3. The nodal value refers to the GQ order applied to either the LWP region when either liquid-only or a mix of cloud phases was present, or IWP when ice-only clouds were present (see Figure 2). As mentioned earlier, 0 were sampled randomly 10 times for each domain. Figure 3 shows errors in upwelling fluxes at TOA for McICA and the first three orders of GQ-ICA relative to ICA benchmarks. McICA is unbiased and so its mean-bias errors are essentially zero. Its random noise, however, can be large: standard deviations in the Tropics reached ∼3 W/m 2 for LW fluxes and ∼25 W/m 2 for both SW and net fluxes. These are consistent with previous studies (e.g. Räisänen et al., 2005). In contrast, the 1-node GQ-ICA displays clear biases, especially for SW fluxes where it overestimates albedo. Evidently, using just the median of W is too often too far up the cumulative distribution. Note that had the mean of W been used, biases would have been even larger (cf. the homogeneous bias as documented in several chapters in Marshak and Davis, 2005). By the 2-node version, however, biases for SW fluxes are all but erased with corresponding standard deviations reduced, at most latitudes, by a factor of ∼4 relative to McICA's. Modest gains were made going to the 3-node version, thus indicating that the additional computation is likely to be unwarranted.
It is interesting to note that the noise for outgoing LW radiation associated with 1-node GQ-ICA exceeded McICA's by a factor of ∼2, while both the 2-and 3-node versions performed approximately the same as McICA. Hence, the main gain made by GQ-ICA is noise reduction for SW transfer. Bringing the fluxes together, GQ-ICA of all orders reduces TOA net flux noise relative to McICA. Figure 4 is like Figure 3 except it is for downwelling fluxes at the surface (SFC). As is well understood, SW results basically mirror those at TOA. But for LW fluxes, all orders of GQ-ICA reduce noise relative to McICA, and this contrasts with TOA results. The reason for this is that downwelling LW fluxes are determined much by low-level cloud distributions, which are dictated much by LWP, and in turn, often sampled better than IWP with double sorting. Like TOA net fluxes, downwelling net fluxes at SFC display ∼3 times less noise than McICA. Figure 5 is like Figure 3 except it is for total absorption by the atmosphere (ATMOS). Like SW fluxes at TOA and SFC, noise associated with GQ-ICA models is substantially smaller than that for McICA. Once again, this demonstrates  Figure 6 shows mean-bias errors (MBEs) for SW fluxes as functions of latitude and 0 for the various models. The large fluctuations shown by McICA are purely due to sampling. Once averaged over 0 they disappear, as seen in Figure 3. Plots for 1-node GQ-ICA, however, echo Figure 3, emphasizing that its biases are independent of 0 . By 2-nodes the biases have been largely eliminated with only slight biases at large 0 . The minor re-emergence of biases at large 0 in the Tropics for 3-node GQ-ICA is because the atmospheres acted on by the GQ-ICA and ICA differed slightly due to double sorting adjusting the location of ice relative to liquid cloud. Placing sub-columns with the largest values of IWP over those with the largest LWPs preserves mean optical depth but increases its variance. This on its own produces domains that are slightly less reflective than ICA. Figure 7 shows root mean square errors (RMSE) that correspond with Figure 6. Throughout the Tropics McICA values at large 0 exceed 30 W/m 2 for TOA and SRF fluxes. Near diurnal-mean 0 they are 25-35 W/m 2 (cf. Figures 3 and 4). It is immediately obvious that GQ-ICA values of RMSE are much smaller than McICA's for all 0 , and that for 2-and 3-node versions they are 10-15 W/m 2 for the most variable tropical clouds at large 0 . This is quite impressive given that representations of cloud-overlap structure get left to just two or three sub-column samples. Figures 7 and 8 show that when it comes to total SW absorption by the atmosphere, GQ-ICA improves on McICA less so than at TOA and SFC. While bias errors are almost negligible everywhere, RMSEs for 2-node GQ-ICA are about 1/2 to 1/3 McICA's values. This might lead one to expect heating rate (HR) RMSEs for 2-and 3-node GQ-ICAs to be substantially smaller than McICA's, but Figure 8 shows that this is not the case. For the SW HRs, RMSEs for 2-node GQ-ICA are only slightly smaller than McICA's, while for the LW those for even 3-node GQ-ICA exceed McICA's. The upshot is that McICA's RMSEs for net HR were the smallest at all altitudes.
The reason for this is that GQ-ICA's very sparse sampling of cloud vertical structure can lead to much uncertainty FIGURE 7 As in Figure 6 except these are root mean square errors in the positioning of cloud tops exposed to space, and this can, in turn, hamper estimation of HR profiles. Hence, while GQ-ICA suitably samples LWP and IWP, there is no constraint on cloud vertical position. Since HR maximizes at cloud tops, getting their positions wrong, even by a single layer, can have a large impact on HR profiles with little impact on corresponding boundary fluxes. With McICA, on the other hand, the many sampled sub-columns yield much better chances at positioning cloud tops well with less uncertainty on the vertical locations of peak HRs. When it comes to HR MBEs, there is little to comment on as all models performed extremely well.
Computational efficiency of radiative transfer algorithms can be important for GCMs. For an atmosphere with a single layer of homogeneous cloud, computation time for 1-node GQ-ICA will be almost equal to McICA's. In general, however, 1-node GQ-ICA's sorting process(es) and the possibility of additional calculations for ice-only sub-columns for double sorting, conspire to make its CPU requirement, for this dataset, ∼1.3 times more than McICA's. For the 2-and 3-node versions the requirements are, respectively, ∼1.8 and ∼2.3 times more than McICA's. Nevertheless, the 2-node GQ-ICA should be comparable to, or faster than, Shonk and Hogan's (2010) tripleclouds scheme for it requires three full-spectrum calculations per cloudy domain.
To be fair to McICA, its CPU usage could be increased to roughly match those of n G -node GQ-ICA by drawing additional sub-column samples for specific CPS points in the CKD model that contribute much to cloud radiative effects. According to fig. 7 in , if regular McICA, as used in this study, were to sample ∼1.8 K, as opposed to K, sub-columns, RMSEs for reflected SW flux at TOA could be reduced from peak diurnal-mean values of ∼25 W/m 2 (see Figures 3 and 6) to ∼10 W/m 2 , which is comparable to those for 2-node GQ-ICA. Nevertheless, McICA's SW and LW routines would still operate on different atmospheric samples. It remains to be seen if this is an issue for large-scale models.

SUMMARY AND CONCLUSION
Many large-scale atmospheric models employ the McICA method for computing domain-average broadband radiative fluxes. McICA has been shown to yield unbiased estimates, relative to the ICA, with accuracies that depend on the quality K/day K/day K/day FIGURE 8 Mean, median, and 0.15 and 0.85 quantiles of heating rate differences between the model listed above each row of plots and ICA benchmarks for all CloudSat-CALIPSO 256 km domains and sampled 0 as described in section 4. Hence, these represent global integrations of the stochastic sub-column cloud generator Pincus et al., 2006;Barker, 2008). It does, however, have the potential to produce substantial amounts of radiative noise at the inner-scale of the host dynamical model. While this noise seldom creates problems for simulations of weather and climate (e.g. Barker et al., 2008), the fact that corresponding SW and LW radiation schemes do not, in general, operate on the same sample atmospheres is ideologically dissatisfying. In this study, a version of McICA that utilizes Gaussian quadrature (GQ) integration -referred to as GQ-ICA -was proposed. GQ-ICA differs from the regular McICA in that the SW and LW models operate on the same cloudy sub-columns.
In essence, GQ-ICA takes several sub-columns, as does McICA, sorts them on cloud water path from smallest to largest, selects sub-columns from the sorted set based on well-defined positions as dictated by GQ, applies full SW and LW spectral integrations to them, weights the results, again from GQ theory, and sums them. Thus, both SW and LW models operate on the same stochastic sample. In GQ integration, the integrand (i.e. fluxes from the transfer models) has to be relatively smooth. For the most part this is true; so much so that even the 2-node GQ-ICA yields almost unbiased flux estimates with substantially less random noise than regular McICA.
While several sorting procedures were tried, most of the results presented were for the so-called double sorting method. This method recognizes that the sum of cloud LWP and IWP cannot be sorted well; LWP is often much greater than IWP, and can easily lead to poor, and even no, GQ sampling of IWP. Of the sorting methods tried, the best performer requires additional computation in that an additional 1-node calculation is often required. For SW radiation, radiative noise for boundary fluxes produced by even 1-node GQ-ICA was much less than McICA's, but it displayed significant mean-bias error, which in turn, was effectively eliminated by the 2-node version. For LW radiation, however, the 2-node GQ-ICA was required to produce noise levels similar to McICA's. Random noise for net radiative fluxes at the surface and TOA was 1/2 to 1/3 that of the regular McICA.
What all GQ-ICA models, for n G ≤ 3, did not do was reduce atmospheric heating rate (HR) noise relative McICA. This is because unlike cloud water path, GQ-ICA lacks direct control over vertical location of sampled cloud. Since GQ-ICA relies on just n G sampled sub-columns, cloud altitude can be sampled poorly, especially when layer cloud fractions are small, thus leading to uncertain HRs. This is most important for cloud tops that are exposed to space. Moreover, it is bothersome given that Räisänen et al. (2005;2008) found HR noise, not surface irradiance noise, to be responsible for very noisy renditions of McICA spawning statistically significant impacts on GCM simulations. At the time of writing, methods to reduce HR noise for GQ-ICA were being investigated. Nevertheless, the high accuracy and computational efficiency of GQ-ICA is encouraging for it suggests that it might be suitable for use in weather and climate models; especially as other taxing sub-models enter the fray, such as circulating oceans and aerosol chemistry, thus reducing the relative amount of CPU resources devoted to radiative transfer.

APPENDIX: ALTERNATIVE SORTING METHODOLOGIES
During the course of this study, several sorting methods, that lead up to Equation 6, were tried. The most successful was the double sorting method that was used to produce results shown in section 5. This appendix explains and shows results for the other sorting methods. Figure A1a is a repeat of Figure 2a. It shows the cells that contained cloud as they came out of the generator. As mentioned in the main text, the simplest way to proceed with sorting when liquid and ice clouds are present in a column is to sort the sub-columnar sums of LWP + IWP. This is referred to as combined sorting and is shown in Figure 2b. First, ∼33% of the sub-columns contained ice but no liquid, ∼5% contained liquid but no ice, and ∼53% contained both phases. Hence, this domain's sub-columns were sorted into two groups as indicated. Computation of F c requires at least two applications of a 1-node GQ; one for each sorted range. By definition, IWP increases monotonically in region-1, but in region-2 it is LWP + IWP that increase monotonically; in general, LWP and IWP alone in region-2 are not expected to increase monotonically. Figure A2a,b show that combined sorting is relatively inferior to double sorting, shown if Figure 2g,h, and what was used for section 5, as there were many large flux errors (random and bias) relative to ICA.
Since high-altitude ice clouds often have large impacts on LW fluxes, another option is to sort only IWP with liquid clouds going along passively. This is referred to as pro-ice cloud sorting and is shown in Figure A1c. While IWP was sorted well, LWP values are almost random. This follows from the two clouds being separated by ∼7 km and the generator treating them almost as random overlap. Note that columns with only liquid clouds were sorted for LWP. As such, for this scenario an n G -node GQ gets applied to region-2 with an additional 1-node GQ to region-1. Figure A2 shows that compared to combined sorting, LW upwelling fluxes at TOA for pro-ice cloud sorting become Errors in 2-node GQ-ICA relative to ICA benchmarks for upwelling SW and LW fluxes at TOA for various LWP and IWP sorting schemes as explained in the text and illustrated in Figure 2. Each point represents results for a randomly sampled 0 using 150 stochastically generated sub-columns based on cloud properties derived from A-Train data. For each domain 10 0 were used. is the root mean square error. Also indicated are percentages of cases with absolute errors >5 W/m 2 slightly more accurate but reflected SW becomes poor with ∼44% of cases showing differences relative to ICA greater than 5 W/m 2 with standard deviation of errors increasing by a factor of 4. This is because LWP was sampled almost randomly, and with just a very small number of samples, two in these cases; this sorting process resembles closely what Räisänen et al.'s (2005) referred to as 1COL-McICA.
Conversely, one can sort LWP with ice cloud going along passively. This is referred to as pro-liquid cloud sorting. Figure A1d shows that outside the liquid cloud region there are sub-columns with ice clouds only. These ice clouds get sorted on IWP. As expected, and as Figure A2 shows, estimates of SW fluxes are slightly better than for combined sorting, but now LW flux estimates are quite poor, but not as poor as SW fluxes were for pro-ice cloud sorting, owing to the majority of IWP values being sampled almost at random. LW transfer depends much on local temperature and is therefore not so affected by random cloud selection as are SW fluxes.
The general implication, therefore, is that sorting of LWP (IWP) should take priority in order to accurately estimate SW (LW) fluxes. Moreover, the double sorting method that was used to produce results shown in section 5 clearly reduces SW and LW flux noise levels simultaneously best.