Assessing CLUBB PDF Closure Assumptions for a Continental Shallow‐to‐Deep Convective Transition Case Over Multiple Spatial Scales

Assumed‐PDF (probability density function) higher‐order turbulence closures (APHOCs) are now widely used for parameterizing boundary layer turbulence and shallow convection in Earth system models (ESMs). A better understanding of the resolution‐dependent behavior of APHOCs is essential for improving the performance of next‐generation ESMs with intended horizontal resolutions finer than 10 km. In this study, we evaluate the PDF family of Analytic double‐Gaussian 1 implemented in Cloud Layers Unified By Binormals (CLUBB) over a range of spatial scales (Dx) from 2 to 100 km. A 120‐km‐wide large eddy simulation (LES) for a continental convection case during 2016 Holistic Interactions of Shallow Clouds, Aerosols, and Land‐Ecosystems (HI‐SCALE) field campaign serves as benchmark to evaluate the PDF closure using an off‐line approach. We find during the shallow convection period, the CLUBB PDF closure tends to produce positive biases of cloud properties and liquid water flux near cloud base for all scales of analysis. It produces negative biases for these variables near cloud top that are more severe for Dx larger than 25 km. Results show that replacing the CLUBB‐parameterized moisture and temperature skewnesses with LES‐derived ones can fix most of the biases if clipping of input moments is allowed to prevent the occurrence of unrealizable solutions. Overall, the performance of the PDF closure is better for smaller Dx = 2–5 km than for larger Dx = 50–100 km; for a given grid spacing, it is better when the convective clouds become deeper in the late afternoon. Likely causes for the resolution dependence and implications for improving the PDF closure are discussed.


Introduction
Boundary layer clouds and turbulence play an important role in regulating the energy and hydrological cycles of the atmosphere but involve processes that are too small in spatial and temporal scales to be explicitly resolved by current or even next generation Earth system models. The conventional approach to representing the cloud-topped boundary layer in Earth system models is to parameterize planetary boundary layer (PBL) processes and moist convection using separate modules. In reality, the processes involved are not always distinct, so the preferred solution is a unified treatment of turbulence and moist convection in atmospheric models (Arakawa, 2004).
The non-Gaussian assumption on the PDF shape (or the PDF closure in short) is what distinguishes APHOCs from traditional third-order closures. It is able to account for the highly skewed variability, which is typical of moist convection. But it is not the only assumption used in closing the APHOC prediction equations. Additional assumptions, similar to those used in traditional second-order closures, are used to parameterize the pressure-correlation and dissipation terms (André et al., 1978;Launder, 1975). APHOCs also make simplifying assumptions to remove terms representing "horizontal interactions" in the prediction of second-and third-order moments (see Shi et al., 2019, for a recent discussion) in their current implementations. All these assumptions are responsible for the overall resolution-dependent behavior of an APHOC scheme.
In this study, instead of looking at the overall behavior of APHOC schemes in a single column or global model (e.g., Guo et al., 2014Guo et al., , 2015Zhang et al., 2018), we isolate the PDF closure from other parts of the APHOC assumptions and evaluate its resolution-dependent behavior using output from a large eddy simulation (LES) with a 120-km-wide domain. We focus on the performance of the Analytical Double Gaussian 1 (ADG1) PDF closure described in Larson and Golaz (2005, referred to as the CLUBB PDF closure hereafter) for a continental convection case with shallow-to-deep convective cloud transitions. The general approach in this study is similar to that in Larson et al. (2002) and Bogenschutz et al. (2010). Note that the early version of the ADG1 PDF family described in Larson et al. (2002) and used in Golaz et al. (2002) is slightly different from what we implement in this study. We would mention the differences in section 3. The early version was tested in Bogenschutz et al. (2010) and implemented in the Simplified Higher-Order Closure (SHOC, Bogenschutz & Krueger, 2013). The detailed analysis using high-frequency output from the large-domain LES corroborates some of the previous conclusions but also exposes problems unnoticed so far. First, when we increase the horizontal scale of analysis beyond 25 km, the overall performance of the CLUBB PDF closure during the shallow convection period deteriorates significantly. Second, we find persistent overestimation of mean cloud properties (cloud fraction, cloud water, and cloud water flux) near cloud base across all horizontal scales we examined and an underestimation of them near cloud top that worsens as the horizontal scale of analysis increases during the shallow convection period. Most of these biases can be fixed by replacing the CLUBB-parameterized moisture and temperature skewnesses with LES-derived ones if clipping of input moments is allowed to prevent the PDF closure from producing unrealizable solutions. This paper is organized as follows. Section 2 introduces the benchmark LES we use and our basic methodology. In section 3 we discuss briefly several assumptions made in the CLUBB PDF closure and introduce other PDF closures we tested alongside the CLUBB PDF closure. Sections 4 and 5 present the results of our analysis that focus on the resolution dependence of CLUBB PDF closure for shallow clouds and its transition to deeper convection respectively. Section 6 concludes and discusses the implications of our study for APHOC application in high-resolution models.

Data and Basic Methodology
The LES we use is based on a daytime continental convection case observed on 30 August 2016 during the Holistic Interactions of Shallow Clouds, Aerosols, and Land-Ecosystems (HI-SCALE) field campaign around the Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) site in north central Oklahoma . The case was simulated using the Weather Research and Forecasting (WRF) model in a LES mode. One-way nesting is used with a 297-km-wide outer domain (Δx = 300 m) and a 120-km-wide inner domain (Δx = 100 m). The initial and time-dependent boundary conditions for the outer domain are obtained from the National Center for Environmental Prediction Final (FNL) operational model global tropospheric analyses (National Center for Environmental Prediction, 2000). The analyses are available at 6-hr intervals on a 1 • × 1 • grid. The vertical resolution is around 24 m below 6 km. The simulation period is from 0600 to 1800 CST (1200-2400 UTC) 30 August. The model three-dimensional output is saved every minute. Setup and detailed evaluation of the WRF-LES simulations can be found in . As a brief summary, it is found that the model can capture the thermodynamical structure and turbulence statistics (e.g., vertical velocity statistics) in the boundary layer and cloud layer reasonably well. When observed high-resolution soil moisture distribution is incorporated into the initial condition, the model can also produce spatial variability in the shallow cloud field that mimics satellite observation. Our analysis is based on the simulation with observed high-resolution soil moisture initial condition, named "revised-1" in .
In this study, we use output only from the inner domain (Δx = 100 m). To minimize the influence of boundary inflow from the outer domain (Δx = 300 m), we exclude grid points that are within 10 km of the lateral boundary, resulting in a 100-km-wide analysis domain. Figure 1 shows the temporal evolution of the domain-averaged cloud fraction and nonprecipitating cloud condensate. Shallow cumulus starts to form in parts of the domain at 1000 CST and transitions to deeper, precipitating convection after 1600 CST. While the primary focus of our study is on the shallow convection period from 1200 to 1400 CST, we also analyze results from 1600 to 1800 CST, which we refer to as the shallow-to-deep transition period. The horizontal distributions of liquid water path averaged over those two periods are depicted in Figure 2. We note that the isolated cloud band appearing after 1500 CST at around 6 km in Figure 1 originates from the outside.
To analyze the resolution-dependent behavior of the CLUBB PDF closure, we first divide our 100-km-wide analysis domain into equal-size subdomains. This is similar to the method used in Jung and Arakawa (2004) (see also Bogenschutz et al., 2010). We denote the subdomain size as Dx = 2, 5, 12.5, 25, 50, and 100 km. Take Dx = 2 km, for example, there are 20 × 20 LES grid points in each 2-km-wide subdomain and in total 50 × 50 subdomains at each vertical level.
Then, the statistical moments used for input to the PDF closure and for evaluation of the output from the PDF closure are calculated for each subdomain as a simple average: where n is the total number of LES grid points within the subdomain and x and y can represent variables like liquid water potential temperature l , total water mixing ratio q t , vertical velocity w, or cloud water mixing ratio q l , and l, m are integers. Cloud fraction (C) is defined as the fraction of cloudy LES grid points (q li > 0) within the subdomain.
Finally, we pass the LES-derived input moments to the PDF closure, diagnose the cloud properties (cloud fraction and cloud water) and the higher-order moments needed in prediction equations using CLUBB PDF closure, and compare those outputs with ones directly computed from LES benchmark for quantitatively analyzing the resolution-dependent behavior of CLUBB PDF closure. In this study, we use our own Python implementation of CLUBB PDF closure and other closures we describe in section 3.

CLUBB PDF Closure
The CLUBB PDF closure uses a double-Gaussian PDF family. The joint PDF, P( l , q t , w), is given by a mixture of two trivariate Gaussians, G 1 and G 2 : where a is the relative weight of the first Gaussian component (G 1 ). One needs 19 parameters (a plus 9 for each Gaussian component) to fully determine the double-Gaussian PDF form, while CLUBB only predicts 10 statistical moments, namely, three first moments (̄l,q t , andw), six second moments , and one third moments (w ′3 ). To solve for the 19 unknowns in terms of the 10 knowns, additional simplifying assumptions are made. We briefly discuss several of these assumptions that are relevant for our analysis. Please refer to Larson and Golaz (2005) for more details.
(i) The CLUBB PDF closure assumes that the variances of w in G 1 and G 2 , 2 w 1 and 2 w 2 , are equal and parameterizes them as̃2 where is a dimensionless constant between 0 and 1.
are the correlation coefficients. The tilde symbol (%) denotes the normalized variable. Note that the early version of the CLUBB PDF closure introduced in Larson et al. (2002) had the variances of w in both Gaussian components set to the same constant (̃2 w 1, 2 = 0.4).

10.1029/2020MS002145
(ii) Earlier work by Lewellen and Yoh (1993) proposed a method to determine the weight of a using the maximum of l , q t , and w skewnesses (denoted by s l , s q t , and s w , respectively, hereafter). The CLUBB PDF closure, however, determines a using only s w (ŝ w to be exact, see below for definition), like in the double delta (DD) function PDF: The correlations between w and q t ( l ) in both G 1 and G 2 are assumed to be 0. (iv) Further, we can derive the following expressions for the mean values of q t in G 1 and G 2 from the We note that assumptions (i) and (iii) are used in the derivation (see Larson & Golaz, 2005, for further details). Here we see that the difference betweenq t1 andq t2 , or the partition of moisture between G 1 and G 2 , does not depend on s q t but onŝ w (in the expression for a). s q t is only involved in determining the variances of q t in G 1 and G 2 , If the difference between s q t and s w , or more precisely the magnitude of (s q t −ĉ 3 wq tŝ w ), becomes too large, we may have to resort to clipping to keep̃2 q t1 and̃2 q t2 nonnegative. Expressions equivalent to Equations 6-9 can be derived for l . Partly to avoid such clipping (Larson & Golaz, 2005, see their discussion in Section 3), the CLUBB PDF closure uses an ansatz to parameterize s q t instead of predicting q ′3 t : where is a tunable parameter between 0 and 3. The same is done for s l . In the early CLUBB version , s l was set to 0 while s q t was set to −1.2s w . In the following analysis, we will compare results from the CLUBB PDF closure with those from several other PDF closures (listed in Table 1) to better understand its behavior. The DD function PDF can be seen as a special form of the general double-Gaussian PDF if all the variances in G 1 and G 2 are 0 and a is fully determined by s w only. The single-Gaussian (SG) PDF is another simpler special form when G 2 = 0 and a = 1 in Equation 3. IPHOC (Cheng & Xu, 2006) uses the ADG2 PDF described in Larson et al. (2002). They also added the prediction equations for ′3 l and q ′3 t (which will be derived from the LES in our case) instead of using the parameterization in Equation 10. In CLUBB_sk, we follow IPHOC and replace the parameterized s q t and s l in the CLUBB PDF closure with those directly derived from the LES. In CLUBB_exp[1-3] we alter the default values for and for sensitivity tests.

Cloud Property Mean Biases
The spatially and temporally averaged profiles of cloud fraction, cloud water mixing ratio, and liquid water flux from the CLUBB PDF closure and LES benchmark for Dx = 2, 5, 25, and 50 km during the shallow convection period are shown in Figures 3a-3c. First, we notice that there is a persistent underestimation of cloud base height of several hundred meters in terms of cloud fraction and cloud water for all Dx by the CLUBB PDF closure. The liquid water flux profiles also show errors consistent with the lower diagnostic cloud base. Cloud fraction, cloud water, and liquid water flux are underestimated by the CLUBB PDF closure in the upper cloud layer. In particular, the errors of cloud water and liquid water flux, compared with cloud fraction, depend more on Dx and become smaller in magnitude for smaller Dx. Both CLUBB cloud water and liquid water flux profiles display a double-peak shape, one near the peak level in the LES profiles and one near the LES cloud base level. This double-peak shape disappears completely only when Dx = 2 km (the blue lines). Also, in the upper cloud layer the mean biases in terms of cloud water and liquid water flux is the smallest for Dx = 2 km.
The biases near cloud top and base are quite insensitive to the choice of and (not shown). The underestimation of cloud base height by the CLUBB PDF closure can also be seen in Figure 1 of Bogenschutz et al. (2010) for the BOMEX case even though less pronounced. On the other hand, we find that the cloud fraction and cloud water profiles from IPHOC shows no underestimation of cloud base height (not shown). This prompts us to look at the difference between the LES-derived moisture and temperature skewnesses, which we feed into IPHOC and those parameterized by CLUBB.
As shown in Figure 4, below the LES cloud base, s q t from the LES (green solid lines) decreases to a negative number very quickly with decreasing altitude while the parameterized s q t from the CLUBB PDF closure (green dash lines) is significantly larger than LES-derived s q t and stays positive for both Dx = 2 and 50 km. Couvreux et al. (2007) found that in a continental dry convective boundary layer the "dry tongues" originated from the entrainment zone is responsible for the negative s q t in their LES. The same process is likely responsible for the negative s q t in the subcloud layer in our case. Above the cloud base, for Dx = 50 km, the magnitude of moisture and temperature skewnesses from the LES increases strongly toward the cloud top (>4 at cloud top), while the parameterized ones reach their maxima of ∼2 in magnitude within the cloud layer. The contrast between the values from LES and CLUBB PDF closure is much smaller for Dx = 2 km. Overall, we see a good correspondence in the vertical between cloud property biases and s q t (s l ) biases. To examine the impact of s q t and s l on the mean cloud properties, we run a modified version of the CLUBB PDF closure, that is, CLUBB_sk, which uses LES-derived s q t and s l instead of parameterized ones. The results are shown in Figures 3d-3f. We immediately see an overall improvement in all three variables, especially in the lower cloud layer. The double-peak shape in the cloud water and cloud water flux profiles for Dx > 2 km disappears but the maximum cloud water flux in the cloud layer is now underestimated. The improvement near the cloud top is not as significant as that in the lower cloud layer for large Dx. We only see a modest increase in cloud water and cloud water flux. We also see frequent occurrence of negativẽ2 q t1,2 and/or̃2 l1,2 near cloud top and base (not shown), which invokes clipping. This indicates that the difference Figure 3. Temporally and spatially averaged profiles of (a) cloud fraction (%), (b) cloud water mixing ratio (q l ), and (c) liquid water flux (w ′ q ′ l ) diagnosed from CLUBB, along with corresponding profiles computed from LES, for the shallow convection regime (1200-1400 CST). For clarity, only the results of Dx = 2, 5, 25, and 50 km are shown. While the averaged 100-m LES cloud fraction and q l profiles have no grid spacing dependence, slight decline in the value of subgrid-scale w ′ q ′ l exhibits between Dx = 2 km and others within the cloud layer. Cloud layer where LES q l > 0.01 g kg −1 is denoted by the blue shading. Panels (d)-(f) are same as panels (a)-(c) but for the CLUBB_sk experiment using LES-derived temperature and moisture skewnesses as input moments into CLUBB PDF closure between LES-derived s q t (s l ) and s w is often too large in these altitudes to produce realizable solutions under the CLUBB PDF closure, as discussed in Larson and Golaz (2005, section 3) and under (iv) in section 3.
The response we see in CLUBB_sk near cloud base can be explained as follows. Given that theĉ wq t in Equations 8 and 9 is positive (i.e., upward moisture flux), when the smaller LES-derived s q t is used near cloud base, the q t variance in the first Gaussian component (G 1 ) would decrease. G 1 is the component with positive mean vertical velocity by definition and in which clouds are more likely to form near cloud base due to larger q t1 . The decrease of q t variance in G 1 near cloud base then leads to cloud initiation at a higher altitude and a decrease in cloud water since q t 1 is below saturation most of the time. The same argument applies to explain the response near cloud top. But for large Dx, the impact of clipping near cloud top is much stronger and damps the response much more significantly than near the cloud base in CLUBB_sk. The clipping method we use in CLUBB_sk is (1) set the negative variance in G 1 or G 2 to 0 then (2) reduce the variance in the other Gaussian so that the total variance stays the same as that from the input (i.e., Equations (13) and (15) in Larson & Golaz, 2005, still hold). However, with this clipping method, the s q t (s l ) produced by the off-line PDF closure calculation in CLUBB_sk differs from that in the LES-derived input (see Equations 8 and 9 or Equation 19 in Larson & Golaz, 2005). We performed another experiment similar to CLUBB_sk except we skipped the second step when clipping. With this alternative clipping method, the q ′2 t ( ′2 l ) produced by the off-line PDF closure calculation differs from that in the input but the difference between the s q t (s l ) produced by the PDF closure and that from the input is reduced. We found that with this change, the impact of clipping near cloud top gets reduced significantly. The cloud water response becomes comparable in magnitude to that near cloud base (not shown).
We argue that the cloud property biases near cloud top and base in the original CLUBB experiment and the frequent clipping in CLUBB_sk show that when the skewness of vertical velocity is significantly different from those of temperature and moisture (especially the latter), the CLUBB PDF closure and other similar PDF closures, like the one used in IPHOC, may be able to capture the vertical velocity variability, that is, the updraft-downdraft structure, but they cannot represent the variability of moisture and temperature in a satisfactory way. In the shallow convection period, this seems more severe for large Dx because of the larger difference between s q t (s l ) and s w near cloud top. The difficulty to capture simultaneously the dynamical and thermodynamical structures of clouds in PDF schemes has been discussed before by Zhu and Zuidema (2009) where they also showed that large differences among s w , s q t , and s l exist in other cases of shallow convection.

Resolution-Dependent Behavior in the Cloud Layer
Next we evaluate the resolution-dependent behavior of the CLUBB PDF closure in the shallow convection layer, defined as vertical levels where the mean cloud water mixing ratio averaged over the shallow convection period from the LES is larger than 0.01 g kg −1 . These vertical levels are indicated by the shading in Figure 3 and more or less encompassed by the first solid contour line in Figure 1b for this time period. To give a sense of the sample size of our analysis, we have 7,440 data points for Dx = 100 km and more than 18 million data points for Dx = 2 km. We also note that our conclusions show little sensitivity to the exact threshold value of cloud water mixing ratio we use to define the cloud layer.
During the shallow convection period, the performance of the CLUBB PDF closure, in terms of correlation with the LES-derived variables, is more or less uniform in time; for a given time, the correlation is typical higher in the middle of the cloud layer and drops off toward cloud top and base. To compare across the Dx range, we construct three metrics for the overall performance for a given variable and a given Dx, namely, correlation coefficient, normalized mean bias, and normalized variance bias. For a particular Dx, the correlation coefficient, R, for a given variable X, is defined as where N is the total number of subdomains within our analysis region between 1200 and 1400 CST,X is the mean over all subdomains and X is the standard deviation. The normalized mean bias (referred to as Figure 7 shows the metrics of w′q ′2 t . The correlation coefficient, R, is very close to 1 for Dx = 2 km and decreases to less than 0.8 for Dx = 100 km. It shows some sensitivity to but almost no sensitivity to . The mean and variance biases show much larger sensitivity to and even though the overall increasing trend with increasing Dx is seen in all but one set of parameter choices. For w ′2 q ′ t in Figure 8, the decreasing trend of R with increasing Dx is even more pronounced (<0.2 for Dx = 100 km). The mean and variance biases also show consistent increase with increasing Dx in all the sensitivity experiments. The liquid water flux w ′ q ′ l is a key component for computing the buoyancy flux in CLUBB, which is the main source of turbulent kinetic energy in the cloud layer. Similar to the other two, it shows a marked drop in R (R < 0.5 for Dx = 100 km) as well as an increase in the mean bias (Figure 9). The variance bias of w ′ q ′ l , however, peaks at Dx = 25 km. In summary, we see a deterioration in the performance of the CLUBB PDF closure as Dx increases in all three higher-order moment categories. We offer a simple interpretation of this trend in section 6. Out of the three categories, the "flux of flux" terms seem to perform the worst in terms of R, especially for Dx > 25 km. The expression for the "flux of flux" terms in the CLUBB PDF closure is identical to that in the DD function  Figure 9. Same as Figure 7 but for the liquid water flux w ′ q ′ l PDF in form, due to the simplifying Assumptions (i)-(iii) listed in section 3. It is our speculation that this form may be too simple to accurately parameterize the flux of flux terms. Finally, it is important to note that the overall resolution-dependent behavior of these high-order moments in CLUBB_sk (not shown) is similar to that in CLUBB even though the variables involving the liquid water perturbation behave quantitatively better and the transport of q t variance for Dx > 25 km performs worse.

Shallow-to-Deep Transition (1600-1800 CST)
In this section, we discuss results from the shallow-to-deep transition period during the last 2 hr of the simulation (1600-1800 CST). As can be seen in Figure 1, this time period sees more rapid vertical growth of the convective cloud populations than the shallow convection period. Precipitation and ice microphysics processes are also more active. In contrast, there is no surface precipitation in the shallow convection period. The mean profiles of cloud fraction, cloud water mixing ratio, and liquid water flux can be found in Figure  10. We focus again on a layer with mean cloud water mixing ratio greater than 0.01 g kg −1 , as indicated by the light blue shading in Figure 10, which is roughly from 2.4 to 4.6 km in altitude. Some features in Figures  3a-3c can also be found in Figure 10, for example, underestimation of cloud base height (even though less pronounced), underestimation of cloud fraction near the top of our analysis layer, and better agreement in general between the LES and the CLUBB PDF closure at smaller Dx than at larger Dx. Figure 11 shows the metrics we defined in section 4.2 for the shaded layer in Figure 10 for cloud fraction (C), w ′ q ′ l , w ′ q ′2 t , and w ′2 q ′ t . The metrics for CLUBB_exp [1][2][3] is also shown to demonstrate the sensitivity to and . First of all, the decreasing trend in R with increasing Dx is much weaker than that in the shallow convection period. The performance in term of R is much better in the transition period for Dx > 25 km (R > 0.75 for all Dx and all variables) and, like in the shallow convection period, shows only limited sensitivity to and . This is true even if we extend the altitude range of our analysis to above 6 km. The performance in terms of the mean and variance biases shows a more significant decreasing trend with increasing Dx albeit Figure 10. Same as Figure 3 but for the profiles averaged over 1600-1800 CST recognized as the shallow-to-deep transition regime. Figure 11. Statistics for C (a-c), w ′ q ′ l (d-f), w ′ q ′2 t (h-j), and w ′2 q ′ t (k-m) as same as Figures 7-9. Results from the transition period are shown in red. For comparison reason, the statistics from prior shallow cumulus are also shown here in blue. more sensitive to and . Compared to the mean and variance biases during the shallow convection period, the performance in the transition period is at least as good.

Concluding Remarks
The objective of this study is to investigate how scale adaptive the CLUBB PDF closure is when provided with input moments sampled from increasingly small subdomains in an LES. We evaluated the CLUBB PDF closure for a continental shallow-to-deep transition case across horizontal scales ranging from 2 to 100 km using a nested WRF-LES simulation with a large 100-km-wide model domain. By passing the "perfect" statistical moments obtained from LES benchmark into assumed-PDF closures, we avoid biases from other parts of the parameterization and from the host model.
Our results show that during the shallow convection period, the CLUBB PDF closure tends to produce positive biases of cloud fraction, cloud water, and liquid water flux near the cloud base for all the subdomain sizes. It also produces negative biases for these variables near the cloud top that are more severe for subdomain sizes larger than 25 km. We find that replacing the CLUBB-parameterized moisture and temperature skewnesses with LES-derived ones fixes most of the biases near cloud base but also leads to unrealizable solutions in the CLUBB PDF closure and requires clipping. The impact of clipping near cloud top for large Dx is so strong that it significantly reduces the response in cloud properties to the changes in moisture and temperature skewnesses.
Within the shallow cloud layer we find that the overall performance of the CLUBB PDF closure deteriorates with increasing subdomain size, especially for subdomain sizes larger than 25 km. This is true for all the output variables of the closure except cloud fraction, which shows little sensitivity to grid spacing. The overall resolution dependence in the shallow-to-deep transition period is the same as that in the shallow convection period. But for subdomain sizes larger than 25 km, the performance in the transition period is considerably better. In Bogenschutz et al. (2010), they showed that for the marine low cloud cases the CLUBB PDF closure does not show much resolution sensitivity for horizontal scales smaller than 25.6 km. This is consistent with our results for the shallow convection period. But for horizontal scales larger than 25 km our results show a clear deterioration in performance during the shallow period. This may be partly due to the mesoscale features present in the shallow cloud field in our case (see . The results (below 5 km in altitude) from the deep convection case (the "Giga-LES" case) in Bogenschutz et al. (2010) show little sensitivity for subdomain sizes ranging from 5 to 100 km, which is also consistent with our results from the transition period.
Overall, the CLUBB PDF closure performs better at 2 or 5 km than at 50 or 100 km in our continental shallow convection case. But, although the biases near cloud top diminish as the resolution increases, those near cloud base are persistent for all horizontal scales. One potential remedy for that is predicting instead of parameterizing third moments of moisture and temperature (see Cheng & Xu, 2006). In the meantime, we need to better understand the difficulty in producing realizable solutions in the CLUBB PDF closure when s q t (s l ) is significantly different from s w . We do not yet fully understand the factors contributing to the resolution dependence we saw in the shallow convection period. Since the joint PDF for a smaller subdomain is likely "simpler" than for a larger subdomain, a double-Gaussian PDF like the one used in CLUBB is more likely to provide a better fit for a smaller subdomain. The large increase in the magnitude of the mean q t ( l ) skewness with subdomain size shown in Figure 4 can be seen as an indication of the increase in the complexity of the PDF. We are currently exploring ways to quantify the complexity of the PDF and its change with horizontal scale using statistical methods from, for example, Firl and Randall (2015). The difference in the performance of the CLUBB PDF closure between the shallow and transition periods may be linked to the difference in the complexity of the PDF or subgrid variability as well. In the deeper convection layer during the transition period, the typical horizontal scale of convective motions becomes larger than that during the shallow convection period. For a given subdomain size, that means the PDF of the variability within the subdomain is likely simpler than that during the shallow convection period. Some evidence of this can be seen in the analysis of Firl and Randall (2015). They found that fitting shallow convection variability in a smaller domain requires more complicated Gaussian mixture models than fitting deep convection variability in a larger domain. Finally, we recognize that extension of our analysis to shallow convective cases with different environmental conditions is needed to better understand the variations in the quantitative resolution dependent behavior of the PDF closure.
As mentioned in section 1, other parameterization assumptions in an APHOC can also impact its resolution-dependent behavior. For example, Guo et al. (2014Guo et al. ( , 2015, in their parameter sensitivity studies, found that the parameters governing the strength of dissipation also have a significant impact on the performance of CLUBB in both single column and global model setups. The interactions between APHOCs and other model components, including the dynamical core and other physical parameterizations, will further complicate the picture in a full three-dimensional atmospheric model. Larson et al. (2012) tested the CLUBB parameterization as a subgrid-scale turbulence parameterization in a cloud resolving model for horizontal resolutions of 2, 4, and 16 km. The results from their continental shallow convection case show an overestimation of cloud water near cloud base and an underestimation near cloud top (see their Figure 16) for all three resolutions, qualitatively similar to our Figures 3a and 3b. This indicates that the biases due to the PDF closure might be significant in a full model. In the future we plan to analyze the relative contributions of different parameterization assumptions in an APHOC quantitatively in a single column model or a full three-dimensional model.

Data Availability Statement
The simulation data used in this study is available through .