Subgrid variations of the cloud water and droplet number concentration over the tropical ocean: satellite observations and implications for warm rain simulations in climate models

One of the challenges in representing warm rain processes in global climate models (GCMs) is related to the representation of the subgrid variability of cloud properties, such as cloud water and cloud droplet number concentration (CDNC), and the effect thereof on individual precipitation processes such as autoconversion. This effect is conventionally treated by multiplying the resolved-scale warm rain process rates by an enhancement factor (Eq ) which is derived from integrating over an assumed subgrid cloud water distribution. The assumed subgrid cloud distribution remains highly uncertain. In this study, we derive the subgrid variations of liquid-phase cloud properties over the tropical ocean using the satellite remote sensing products from Moderate Resolution Imaging Spectroradiometer (MODIS) and investigate the corresponding enhancement factors for the GCM parameterization of autoconversion rate. We find that the conventional approach of using only subgrid variability of cloud water is insufficient and that the subgrid variability of CDNC, as well as the correlation between the two, is also important for correctly simulating the autoconversion process in GCMs. Using the MODIS data which have near-global data coverage, we find that Eq shows a strong dependence on cloud regimes due to the fact that the subgrid variability of cloud water and CDNC is regime dependent. Our analysis shows a significant increase of Eq from the stratocumulus (Sc) to cumulus (Cu) regions. Furthermore, the enhancement factor EN due to the subgrid variation of CDNC is derived from satellite observation for the first time, and results reveal several regions downwind of biomass burning aerosols (e.g., Gulf of Guinea, east coast of South Africa), air pollution (i.e., East China Sea), and active volcanos (e.g., Kilauea, Hawaii, and Ambae, Vanuatu), where the EN is comparable to or even larger than Eq , suggesting an important role of aerosol in influencing the EN . MODIS observations suggest that the subgrid variations of cloud liquid water path (LWP) and CDNC are generally positively correlated. As a result, the combined enhancement factor, including the effect of LWP and CDNC correlation, is significantly smaller than the simple product of Eq ·EN . Given the importance of warm rain processes in understanding the Earth’s system dynamics and water cycle, we conclude that more observational studies are needed to provide a better constraint on the warm rain processes in GCMs. Published by Copernicus Publications on behalf of the European Geosciences Union. 1078 Z. Zhang et al.: Subgrid cloud property variations and covariations from satellite

Abstract. One of the challenges in representing warm rain processes in global climate models (GCMs) is related to the representation of the subgrid variability of cloud properties, such as cloud water and cloud droplet number concentration (CDNC), and the effect thereof on individual precipitation processes such as autoconversion. This effect is conventionally treated by multiplying the resolved-scale warm rain process rates by an enhancement factor (E q ) which is derived from integrating over an assumed subgrid cloud water distribution. The assumed subgrid cloud distribution remains highly uncertain. In this study, we derive the subgrid variations of liquid-phase cloud properties over the tropical ocean using the satellite remote sensing products from Moderate Resolution Imaging Spectroradiometer (MODIS) and investigate the corresponding enhancement factors for the GCM parameterization of autoconversion rate. We find that the conventional approach of using only subgrid variability of cloud water is insufficient and that the subgrid variability of CDNC, as well as the correlation between the two, is also important for correctly simulating the autoconversion process in GCMs. Using the MODIS data which have near-global data coverage, we find that E q shows a strong dependence on cloud regimes due to the fact that the subgrid variability of cloud water and CDNC is regime dependent. Our analysis shows a significant increase of E q from the stratocumulus (Sc) to cumulus (Cu) regions. Furthermore, the enhancement factor E N due to the subgrid variation of CDNC is derived from satellite observation for the first time, and results reveal several regions downwind of biomass burning aerosols (e.g., Gulf of Guinea, east coast of South Africa), air pollution (i.e., East China Sea), and active volcanos (e.g., Kilauea, Hawaii, and Ambae, Vanuatu), where the E N is comparable to or even larger than E q , suggesting an important role of aerosol in influencing the E N . MODIS observations suggest that the subgrid variations of cloud liquid water path (LWP) and CDNC are generally positively correlated. As a result, the combined enhancement factor, including the effect of LWP and CDNC correlation, is significantly smaller than the simple product of E q · E N . Given the importance of warm rain processes in understanding the Earth's system dynamics and water cycle, we conclude that more observational studies are needed to provide a better constraint on the warm rain processes in GCMs.

Introduction
Marine boundary layer (MBL) clouds are a strong modulator of Earth's radiative energy budget (Klein and Hartmann, 1993;Trenberth et al., 2009). They can interact with other components of the climate system, such as aerosols and precipitations, in various ways. The feedback of MBL clouds to climate change remains one of the largest uncertainties in our understanding of the climate sensitivity (Bony and Dufresne, 2005;Soden and Held, 2006). Despite their importance in the climate system, simulating MBL clouds in general circulation models (GCMs) has proved to be extremely challenging. A main difficulty is rooted in the fact the typical grid size of GCMs (∼ 100 km) is much larger than the spatial scale of many cloud microphysical processes, and as a result these subgrid scale processes, as well as the subgrid cloud variations, have to be highly simplified and then parameterized as functions of resolved, grid-level variables.
Of particular interest in this study is the warm rain processes in MBL clouds, which have fundamental impacts on the cloud water budget and lifetime. Although in reality it is highly complicated and involves multiple factors, warm rain formation in GCMs is usually parameterized as simple functions of only key cloud parameters. For example, the drizzle in MBL clouds is initialized by the so-called autoconversion process in which the collision-coalescence of cloud droplets gives birth to large drizzle drops (Pruppacher and Klett, 1997). In GCMs, for the sake of efficiency, this process is usually parameterized as a power function of liquid water content (LWC or symbol q c ) and cloud droplet number concentration (CDNC or symbol N c ). One of the most widely used parameterization schemes is developed by Khairoutdinov and Kogan (2000) (KK2000 hereafter), which has the form where ∂q r ∂t is the rain water tendency due to the autoconversion process, q c has the unit of kg kg −1 , and N c the unit of cm −3 . The three parameters C = 1350, β q = 2.47, and β N = −1.79 are derived through a simple least-square fitting of the autoconversion rate results from a large-eddy simulation with bin microphysics that can simulate the processlevel physics. Even though this is highly simplified, the parameterization scheme still faces a great challenge. The calculation of grid-mean autoconversion efficiency requires the knowledge of subgrid distributions of LWC and CDNC, but in the GCMs only grid-mean quantities q c and N c are known and available for use in the computation of the autoconversion rate. As pointed out by Pincus and Klein (2000), for a process f (x) such as autoconversion that is nonlinearly dependent on subgrid variables, x, the grid-mean value f (x) is not equal to the value estimated based on the gridmean x , i.e., f (x) = f ( x ). Mathematically, if f (x) is convex, then f ( x ) < f (x) Larson et al., 2001). To take this effect into account, a parameter E is often introduced in the GCM as part of the parameterization such that f (x) = E · f ( x ). It is referred to as the enhancement factor in many studies as well as this study because E > 1 for a convex function. Such a nonlinear effect is not just limited to the autoconversion process. Some other examples are the plane-parallel albedo bias (Barker, 1996;Cahalan et al., 1994;Oreopoulos and Davies, 1998a), subgrid cloud droplet activation (Morales and Nenes, 2010), and accretion (Boutle et al., 2014;Lebsock et al., 2013).
The value of E is determined primarily by two factors: the nonlinearity of f (x) and the subgrid probability density function (PDF) P (x). Given the same subgrid variation of LWC, i.e., P (q c ), the nonlinear effect impacts the autoconversion parameterization more than it does the accretion because the former is a more nonlinear function of q c than the latter. For the same f (x), a grid box with a narrow and symmetric P (x) would require a smaller E than another grid box with a broader and nonsymmetric P (x). Ideally, the value of the enhancement factor E should be diagnosed from the subgrid cloud PDF P (x). Unfortunately, because this is not possible in most conventional GCMs, the value of E is usually assumed to be a constant for the lack of better options. The E for autoconversion due to subgrid LWC variation is assumed to be 3.2 in the two-moment cloud microphysics parameterization schemes by Morrison and Gettelman (2008) (MG scheme hereafter), which is employed in the widely used Community Atmosphere Model (CAM). This choice of E = 3.2 is based on an early study by Barker et al. (1996), in which the mesoscale variation of column-integrated optical thickness of the "overcast stratocumulus", "broken stratocumulus", and "scattered stratocumulus" are studied. The value E = 3.2 is derived based on the mesoscale variation of the broken stratocumulus.
Clearly, a simple constant E is not adequate. The following is a list of attempts to better understand the subgrid cloud variations and the implications for warm rain simulations in GCMs. Several previous studies have shown that the mesoscale cloud water variation is a strong function of cloud regime -the subgrid cloud water variation of Sc clouds is much different from that of Cu clouds Lee et al., 2010;Oreopoulos and Cahalan, 2005;Wood and Hartmann, 2006). As the first part of a two-part study, Larson and Griffin (2013) first laid out a systematic theoretical basis for understanding the effects of subgrid cloud property variations on simulating various nonlinear processes in GCMs, including not only the autoconversion but also the accretion, condensation, evaporation, and sedimentation processes. In the second part, using cloud fields from a largeeddy simulation (LES), Griffin and Larson (2013) showed that inclusion of the enhancement factor indeed leads to more rainwater at surface in single-column simulations and makes them agree better with high-resolution large-eddy simulations. Recently, using the ground-based observations from three Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) sites, Xie and Zhang (2015) developed a scale-aware parameterization scheme for GCMs to account for subgrid cloud water variation. Also using ARM measurement, Ahlgrimm and Forbes (2016) analyzed the dependence of cloud water variability on cloud regime. Although these previous studies have shed important light on subgrid cloud variation and the implications for GCM, they lack a global perspective because they are only based on limited data (e.g., LES cases, in situ and ground-based measurement). Currently, satellite remote sensing observation is the only way to achieve a global perspective. Using the observations from the spaceborne radar CloudSat, Lebsock et al. (2013) showed that the subgrid cloud water variance is smaller over the Sc region than over the Cu region, and as a result the enhancement factor shows an increasing trend from Sc to Cu regions. They also highlighted importance of considering the subgrid covariability of cloud water and rain water in the computation of the accretion rate. Using a combination of in situ measurement and satellite remote sensing data, Boutle et al. (2014) analyzed the spatial variation of both cloud and rain water, as well as their covariation, and developed a simple parameterization scheme to relate the subgrid cloud water variance to the grid-mean cloud fraction. Later, the study of Boutle et al. (2014) was extended by Hill et al. (2015), who developed a cloud-regime-dependent and scale-aware parameterization scheme for simulating subgrid cloud water variation. On the modeling side, Guo et al. (2014) investigated the sensitivity of cloud simulations in the Geophysical Fluid Dynamics Laboratory (GFDL) atmospheric general circulation model (AM) to the subgrid cloud water parameterization schemes. A similar study was carried out by Bogenschutz et al. (2013) using the National Center for Atmospheric Research (NCAR) CAM. Both studies show that the more sophisticated subgrid parameterization scheme -Cloud Layers Unified by Binormals (CLUBB) (Golaz et al., 2002a, b;Larson et al., 2002) -lead to a better simulation of clouds in the model. However, a more recent study by Song et al. (2018b) reveals that the CLUBB in CAM version 5.3 (CAM5.3) overestimates the enhancement factor in the trade wind cumulus cloud region, which in turn leads to excessive drizzle in the model and "empty clouds" with nearzero cloud water. In addition to CLUBB, the so-called superparameterization (a.k.a. multiscale modeling framework -MMF), which uses cloud resolving models embedded in the GCM grids to diagnose subgrid cloud variations (Randall et al., 2003), has also gained increasing popularity. Takahashi et al. (2017) compared the subgrid cloud water variations simulated by a CAM-MMF model with those derived from A-Train observations and found reasonable agreement.
Despite these previous studies, many questions remain unanswered. First of all, all the previous studies, as far as we know, have focused on the impact of subgrid cloud water q c variation. The potential impact of subgrid variation of N c and the covariability of N c with q c have been overlooked so far. Given the same amount of q c , a cloud with a smaller N c would have larger droplets and therefore larger precipitation efficiency than another cloud with a larger N c . For the same reason, other things being equal, a grid with a positive correlation of subgrid N c and q c would be less efficient in terms of autoconversion than a grid with a negative correlation of the two. Secondly, most of previous studies are based on the assumption that the subgrid cloud property variation follows certain well-behaved distributions, usually either gamma (e.g., Barker, 1996;Morrison and Gettelman, 2008;Oreopoulos and Barker, 1999;Oreopoulos and Cahalan, 2005) or lognormal (e.g., Boutle et al., 2014;Larson and Griffin, 2013;Lebsock et al., 2013). However, the validity and performance of the assumed PDF shape are seldom checked. Furthermore, although the study by Lebsock et al. (2013) has depicted a global picture of the enhancement factor for the autoconversion modeling in GCM, the picture is far from clear due to the small sampling rate of CloudSat observations.
In this study, we revisit the subgrid variations of liquidphase cloud properties over the tropical ocean using 10 years of MODIS cloud observations, with the overarching goal to better understand the potential impacts of subgrid cloud variations on the warm rain processes in the conventional GCMs. Similar to previous studies, we will quantify the subgrid cloud water variations based on MODIS observations. Going one step further, we will also attempt to unveil for the first time the subgrid N c variation, as well as its correlation with cloud water, and investigate the implications for warm rain simulations in GCM. Moreover, we will take advantage of the wide spatial coverage of MODIS data to achieve a more detailed picture of the enhancement factor for the autoconversion simulation. Last but not least, we will evaluate the two widely used distributions, i.e., lognormal and gamma, in terms of their performance and limitations for simulating the enhancement factor. We will first explain the theoretical background in Sect. 2 and introduce the data and methodology in Sect. 3. The MODIS observations will be presented and discussed in Sect. 4. The implications for the autoconversion parameterization in the GCMs will be discussed in Sect. 5. The main findings will be summarized in Sect. 6 with an outlook for future studies. In previous studies, the spatial variations of cloud properties, such as cloud optical thickness (COT), cloud liquid water path (LWP) and cloud LWC, are often described using either of two theoretical distributions -the gamma and lognormal distribution. The PDF from a gamma distribution is a two-parameter function as follows (Barker, 1996;Oreopoulos and Davies, 1998b): where is the gamma function, v is the so-called inverse relative variance, and α the so-called rate parameter. If x follows the gamma distribution, its mean value is given by and variance given by It follows from Eqs. (3) and (4) that the so-called inverse relative variance is where η = Var(x) x 2 is the relative variance. If x follows the gamma distribution for a physical process M (x) that is a power function of x, then the expected value M(x) is given by As explained in the introduction, for a nonlinear process M (x), M(x) = M ( x ). The ratio between the two E values is by definition the enhancement factor: The PDF of a lognormal distribution is given as follows Lebsock et al., 2013): where µ = ln x and σ 2 = Var(ln x) correspond to the mean and variance of ln x, respectively. The mean value of the lognormal distribution is given by and the variance by It follows from Eqs. (10) and (11) that the inverse relative variance can be derived from the following equation: If x follows the lognormal distribution, the expected value of Evidently, the corresponding enhancement factor is given by Note that Eqs. (7) and (8) are only valid when β > −v because gamma function (v + β) can run into singular values when v + β < 0. In contrast, Eqs. (13) and (14) are valid for any real value β. This is one advantage of the lognormal distribution over the gamma distribution. An example of the gamma and lognormal distributions for q c is shown in Fig. 1a. In this example, both distributions have the same mean q c = 0.5 g kg −1 and also the same inverse relative variance v q = 3. Although the general shapes of the two PDFs are similar, they differ significantly at the two ends; the gamma PDF is larger than lognormal PDF over the small values of q c , and the opposite is true over the large values of q c . The gamma and lognormal distributions can also be used to describe the spatial variation of N c ( Gultepe and Isaac, 2004). An example is given in Fig. 1c, in which q c is a constant of 0.5 g kg −1 , N c = 50 cm −3 , and v N = 5.0. Figure 1b shows the autoconversion rate based on the KK2000 parameterization scheme for the gamma P G (q c ) and lognormal P L (q c ) that are shown in Fig. 1a. Interestingly, although the cumulative autoconversion rates based on the two types of PDFs are almost identical, the contributions to the total autoconversion rate from the different LWC bins are quite different. As shown in Fig. 1a, the P L (q c ) has a longer tail than the P G (q c ); i.e., the occurrence probability of large q c (e.g., q c > 2.0 g kg −1 ) is much higher in the lognormal than in gamma PDF. This difference is further amplified in the autoconversion rate computation in Fig. 1b because the autoconversion rate is proportional to q 2.47 c . The enhancement factors based on the gamma (i.e., E (P G , β) in Eq. 8) and lognormal (i.e., E(P L , β) in Eq. 14) PDFs for β q = 2.47 are plotted as a function of the inverse relative variance v in Fig. 2. When subgrid clouds are more  homogenous, e.g., v > 1, the enhancement factor based on the two PDFs are similar. However, for more inhomogeneous grids, e.g., v < 1, the E(P L , β) is significantly larger than E(P G , β), which is probably because of the longer tail of P L (q c ) as shown in Fig. 1a and b.

Impacts of subgrid cloud variations on warm rain parameterization in GCM
The warm rain process in MBL clouds involves many interacting microphysical processes. In this study, we focus only on the simulation of autoconversion in GCM. Other nonlinear processes, such as accretion and evaporation, have been investigated in previous studies (Boutle et al., 2014;Lebsock et al., 2013). Ideally, if the subgrid variations of q c and N c are known, then the grid-mean in-cloud autoconversion rate should be derived from the following integral: where P (q c , N c ) is the joint PDF of q c and N c . Unfortunately, most conventional GCMs lack the capability of predicting the subgrid variations of cloud properties, with only a couple of exceptions (Thayer-Calder et al., 2015). What is known from the GCM is usually the in-cloud grid-mean values q c and N c . As a result, instead of using Eq. (15), the autoconversion rate in GCMs is usually computed from the where E is the enhancement factor defined as The value of the enhancement factor depends on the subgrid variations of q c and N c . If clouds are homogenous on the subgrid scale, then E ∼ 1. In the special case where q c and N c are independent, then the joint PDF P (q c , N c ) becomes P (q c , N c ) = P (q c )P (N c ), where P (q c ) and P (N c ) are the PDF of the subgrid q c and N c . Consequently, Eq. (15) reduces to and Eq. (17) to where E q is the enhancement factor due to the subgrid variation of cloud water, which has the form and the E N is the enhancement factor due to the subgrid variation of CDNC, which has the form Obviously, if P (q c ) and P (N c ) follow either a gamma or lognormal distribution, then the above equations reduce to Eqs. (8) or (14), respectively. If q c and N c both have significant subgrid variations and they are not independent, the enhancement factor should ideally be diagnosed from Eq. (17). However, the joint PDF P (q c , N c ) may not be known, and the integration can be time-consuming. Some previous studies proposed to approximate the P (q c , N c ) as a bivariate lognormal distribution as follows: where ρ is the correlation coefficient between q c and N c Lebsock et al., 2013). As such, both q c and N c follow a marginal lognormal distribution in Eq. (9). Substituting Eq. (22) into Eq. (17), we obtain the enhancement factor for the bivariate lognormal distribution that consists of three terms: where correspond to the impacts of subgrid q c and N c variance, respectively (i.e., Eq. 14), and corresponds to the impact of the covariation of q c and N c on the enhancement factor. Obviously, Eq. (23) reduces to Eq. (19) when q c and N c are uncorrelated (i.e., ρ = 0, E COV = 1). If q c and N c are negatively correlated (i.e., ρ < 0 and E COV > 1), clouds with a larger q c would tend to have a smaller N c . The autoconversion rate in such a case would be larger than that in the case where q c and N c are positively correlated (i.e., ρ > 0 and E COV < 1). A positive correlation would exist, for instance, if all droplets in clouds were the same size, but some parcels had more droplets than other parcels. Most current GCMs do not have the capability to simulate the subgrid cloud property variations. They usually have to use predefined subgrid cloud variations in the computation of grid-mean autoconversion rate instead of using prognostic values. For example, in the MG scheme for the CAM5.3, the subgrid q c is assumed to follow the gamma distribution in Eq. (2) with a fixed v q = 1 and, as a result, a constant E q = 3.2. Lately, advanced subgrid parameterization schemes, such as CLUBB, have been implemented in several GCMs, including CAM6 and GFDL AM model (Bogenschutz et al., 2018;Guo et al., 2015Guo et al., , 2014, which provides information on the subgrid q c variation to the host model. The information can then be used to dynamically diagnose the enhancement factor E q , which will help the model simulate the cloud regime dependence of E q (Guo et al., 2010(Guo et al., , 2014. However, as explained above, not only the subgrid variation of q c but the subgrid variation of N c can also influence the enhancement factor. Unfortunately, this aspect has been ignored by almost all GCMs, even the latest CAM6 with CLUBB. Physically, provided the same q c , a cloud with a smaller N c would have a larger droplet size and therefore larger precipitation efficiency than the cloud with a larger N c . Because the autoconversion rate depends nonlinearly on N c , the grid-mean autoconversion rate computed based on a skewed PDF of N c (i.e., ∞ 0 (N c ) β N P (N c ) dN c ) would be different from that computed based on the mean of N c (i.e., ( N c ) β N ). The autoconversion enhancement factor based on the lognormal PDF E(P L , β) for β N = −1.79 is given in Fig. 2. Interestingly, at the same inverse relative variance v, the enhancement factor based on the same lognormal PDF E(P L β) for β N = −1.79 is actually larger than that for β q = 2.47 because of the formula of the exponent in Eq. (14) (i.e., β 2 −β 2 ). Moreover, the correlation between N c and q c can also be important. Going back to Eq. (23), evidently, E > E q if and only if E N · E COV > 1. After some manipulation, we can show that if β N < 0 and σ N > 0, then This equation reveals that when q c and N c are weakly or negatively correlated (ρ ≤ 0), considering only E q would tend to underestimate E. On the other hand, however, if q c and N c are highly positively correlated (ρ ∼ 1), then considering E q only would tend to overestimate E.

Data and methodology
To derive the above-mentioned enhancement factors, we will use 10 years (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) of the latest collection 6 (C6) daily mean level-3 cloud retrieval product from the Aqua MODIS instrument (product name MYD08_D3), which contains the gridded statistics of cloud properties computed from pixellevel (i.e., level-2) retrievals. As summarized in Platnick et al. (2003Platnick et al. ( , 2017, the operational level-2 MODIS cloud product provides cloud masking (Ackerman et al., 1998), cloud top height (Menzel et al., 1983), cloud top thermodynamic phase determination (Menzel et al., 2006), and COT, cloud effective radius (CER), and LWP retrievals based on the bispectral solar reflectance method (Nakajima and King, 1990). All MODIS level-2 atmosphere products, including the cloud, aerosol, and water vapor products, are aggregated to 1 • × 1 • spatial resolution on a daily, 8-day, and monthly basis. Aggregations include a variety of scalar statistical information, including mean, standard deviation, and max and min occurrences, as well as histograms including both marginal and joint histograms. For COT, CER, and LWP, the MODIS level-3 product provides both their in-cloud gridmean values ( x ) and subgrid standard deviations (σ x ). The inverse relative variance v can then be derived from Eq. (5), i.e., v = x 2 /σ 2 x . Note that the operational MODIS product provides two CER retrievals: one based on the observation from the band 7 centered around 2.1 µm and the other from band 20 at 3.7 µm. As discussed in several previous studies (Cho et al., 2015;Zhang and Platnick, 2011;Zhang et al., 2012Zhang et al., , 2016, the 3.7 µm band CER retrieval is more resilient to the 3-D effects and retrieval failure than the 2.1 µm band retrievals. For these reasons, it is used as the observational reference in this study.
Given the COT and CER retrieval, the operational MODIS product estimates the LWP of clouds using where ρ w is the density of water. Several studies have argued that a smaller coefficient of 5/9, instead of 2/3, should be used in estimation of LWP (Lebsock et al., 2011;Seethala and Horváth, 2010;Wood and Hartmann, 2006). The choice of coefficient does not matter in this study because it is a common factor in the calculation of v. The choice of the coefficient has no impact on our study because we are interested in the relative inverse variance v = x 2 /σ 2 x . We note here that it is the LWC q c , instead of the LWP, that is used in the KK2000 scheme. So, the spatial variability of LWC is what is most relevant. However, the remote sensing of cloud water vertical profiles from satellite sensors for liquid-phase clouds is extremely challenging even with active sensors. It is why most previous studies using the satellite observations analyzed the spatial variation of LWP rather than LWC. In fact, even Lebsock et al. (2013), who used the level-2 Cloud-Sat observations, had to use the vertically averaged LWC in their analysis. Airborne in situ measurement faces a similar challenge. For example, Boutle et al. (2014) used the LWC observation along horizontal flight tracks to study the spatial variability of cloud water, which only samples the LWC at certain levels of MBL clouds. Ground-based observations are much better than satellite and airborne observations in this regard. Recently, Xie and Zhang (2015) analyzed the cloud water profiles retrieved using ground-based radars from the three ARM sites and found no obvious in-cloud vertical dependence of the spatial variability of LWC. Following these previous studies, we assume that the horizontal subgrid variation of LWC is not strongly dependent on height, and its value can be inferred from the spatial variability of the vertically integrated quantities of LWP. The uncertainty caused by this assumption will be assessed in future studies.
The current MODIS level-3 cloud product does not provide CDNC retrievals. Following previous studies (Bennartz, 2007;Bennartz and Rausch, 2017;Grosvenor and Wood, 2014;McCoy et al., 2018), we estimate N c of liquid-phase clouds from the MODIS-retrieved COT (τ ) and CER (r e ) based on the classic adiabatic cloud model where ρ w is the density of water, Q e ≈ 2 is the extinction efficiency of cloud droplets, k is the ratio of r e to mean volumeequivalent radius, f ad is the adiabaticity of the cloud, and w is the LWC lapse rate. Following previous studies, we assume k = 0.8 and f ad = 1.0 to be constant and compute w from the grid-mean liquid cloud top temperature and pressure. The theoretical basis and main uncertainty sources of the CDNC estimation based on the adiabatic cloud model from MODIS-like passive cloud retrievals are nicely reviewed by Grosvenor et al. (2018). Ideally, the values of LWP and CDNC should be estimated on a pixel-by-pixel basis from the level-2 MODIS product. However, pixel-by-pixel estimation is highly timeconsuming, which makes it difficult to achieve a global perspective. Using an alternative method, many previous studies estimate the grid-level CDNC statistics from the joint histogram of COT vs. CER provided in the level-3 MODIS cloud products (Bennartz, 2007;McCoy et al., 2017McCoy et al., , 2018. For a given 1 • × 1 • grid box, the liquid-phase COT-CER joint histogram provides the counts of successful cloud property retrievals with respect to 108 joint COT-CER bins that are bounded by 13 COT bin boundaries, ranging from 0 to 150, and 10 CER bin boundaries, ranging from 4 to 30 µm. With the joint histogram, which is essentially the joint PDF of COT and CER P (τ r e ), we can estimate the grid mean and variance of CDNC from the following equations: x = x (τ, r e ) P (τ, r e ) dτ dr e , where x can be either LWP or CDNC. Figure 3a shows the LWP in Eq. (26) as a function of the 13 COT bins and 10 CER bins from the MODIS level-3 product. As expected, the largest LWP values are found when both COT and CER are large. Figure 3b shows the CDNC in Eq. (27) as a function of the COT and CER bins. As expected, the largest CDNC values are found when both COT is large and CER is small. Figure 3c shows an example of the COT-CER joint histogram from the Aqua MODIS daily level-3 product MYD08_D3 on 9 January 2007 at the grid box 1 • S and 1 • W. In this particular grid box, a combination of ∼ 2-4 COT and ∼ 10-12 µm CER is the most frequently observed cloud value. Using the joint histogram in Fig. 3c, we can derive the mean and variance of both LWP and COT using Eqs. (28) and (29). The efficiency of using the level-3 MODIS product is accompanied by three important limitations. First of all, as mentioned earlier, MODIS provides only LWP retrievals, while LWC is needed in the KK2000 scheme. Second, the current level-3 MODIS cloud product has a fixed 1 • × 1 • spatial resolution. Although this resolution is highly relevant to the current generation of GCMs, i.e., Coupled Model Intercomparison Project Phase 6 (CMIP6) (Eyring et al., 2016), future GCMs may have significantly finer resolution. Third, it is difficult to subsample the pixels with the best retrieval quality. These limitations will have to be addressed in future studies.

Grid-mean and subgrid variations of liquid-phase cloud properties
In this study, we limit our analysis to tropical oceans only where warm rain is frequent and MODIS cloud retrievals have a relatively better quality than over land or over high latitudes. The annual mean total cloud fraction (f tot ), liquidphase cloud fraction (f liq ), in-cloud COT, CER from the 3.7 µm band, LWP, and estimated CDNC over the tropical oceans based on 10 years Aqua MODIS retrievals are shown in Fig. 4. The highest f liq in the tropics is usually found in the stratocumulus (Sc) decks over the eastern boundary of the ocean, e.g., the SE Pacific off the coast of Peru, the NE Pacific off the coast of California, and the SE Atlantic off the coast of Namibia. The liquid-cloud fraction reduces significantly toward the open ocean trade wind regions, where the dominant cloud types are broken cumulus (Cu). Close to the continents, the Sc decks are susceptible to the influence of the continental air mass with a higher loading of aerosols in comparison with a pristine ocean environment, which is probably the reason the SC decks have a smaller CER and higher CDNC than the open-ocean trade cumulus ( Fig. 4d and f). The in-cloud COT (Fig. 4c) and LWP (Fig. 4e) generally increase from the Sc decks to the open-ocean Cu regime, although less dramatically than the transition of cloud fraction. The Sc decks and the Sc-to-Cu transition are the most promi- nent features of liquid-phase clouds in the tropics. However, as mentioned in the introduction, simulating these features in the GCMs proves to be an extremely challenging task, and most GCMs suffer from some common problems, such as the "too few, too bright" problem and the abrupt Sc-to-Cu transition problem (Kubar et al., 2014;Nam et al., 2012;Song et al., 2018a). Switching the focus now from grid-mean values to subgrid variability, we will show the grid-level inverse relative variances v = x 2 /Var(x)) for several key cloud properties. Here, we first derive the daily mean v and then aggregate the result to monthly mean values. Therefore, for each grid box we have 120 samples (i.e., 10 years × 12 months) of monthly mean v for analysis and visualization. Because the value of v can be ill behaved when Var (x) approaches zero, instead of the mean value, we plot the median value of v based on 120 months of MODIS observations in Fig. 5. There are several interesting and important features in Fig. 5. First of all, the v values of all four sets of cloud properties (i.e., COT, CER, LWP, and CDNC) all exhibit a clear and similar Scto-Cu transition, with larger values in the Sc regions and smaller values in the broken Cu regions. This indicates that cloud properties, including both optical and microphysical properties, are more homogenous, in terms of spatial distribution within the grid, in the Sc region than in the Cu region.
Secondly, the value of v of CER (i.e., 10-100 in Fig. 5b) is larger than that of the other properties (i.e., 1-10) by almost an order of magnitude, indicating that the subgrid variability of CER is very small. On the other hand, however, it is important to note that the v of CDNC (Fig. 5d) is comparable with that of COT (Fig. 5a) and LWP (Fig. 5c). The reason is probably in part because the highly nonlinear relationship between CDNC and CER (i.e., N c ∼ r − 5 2 e ) leads to a stronger variability of CDNC than CER and also in part because the variability of CDNC is also contributed by the subgrid variation of COT. In some regions, the Gulf of Guinea, the East and South China Sea, and Bay of Bengal for example, the v of CDNC is close to unity, indicating the subgrid standard deviation of CDNC is comparable to the grid-mean values in these regions. As discussed in the next section, the significant subgrid variability of CDNC in these regions should be taken into account when modeling the nonlinear processes, such as the autoconversion, in GCM to avoid systematic biases due to the nonlinearity effect.
The values of v in Fig. 5 from this study are in reasonable agreement with previous studies. Barker (1996) selected a few dozens of cloud scenes, each about 100-200 km in size, from the Landsat observation and analyzed their spatial variability of COT. It is found that the typical value of v for overcast stratocumulus, broken stratocumulus, and scattered cumulus is 7.9, 1.2, and 0.7, respectively (see their  Table 3), which is consistent with the Sc-to-Cu transition pattern seen in Fig. 5. Oreopoulos and Cahalan (2005) derived the subgrid inhomogeneity of COT on a global scale from the level-3 Terra MODIS retrievals. Although using a different metric (i.e., their inhomogeneity parameter is defined as χ = exp (ln τ ) / τ ), they also found a systematic increase of inhomogeneity (decreasing value of χ) from the Sc region to Cu region. Also using the MODIS cloud property retrievals, Wood and Hartmann (2006) investigated the mesoscale spatial variability of LWP in the NE Pacific and SE Pacific regions. The v of LWP is found to increase systematically with mesoscale cloud fraction, and the relationship between the two can be reasonably explained by a simple PDF cloud thickness model in Considine et al. (1997). See also Kawai and Teixeira (2010).
As explained in Sect. 2, the correlation between cloud water and CDNC can also influence the computation of enhancement factor and thereby the grid-mean autoconversion rate. Figure 5e shows the median value of the LWP and CDNC correlation coefficient ρ. Similar to the derivation of median v, we first compute the monthly mean ρ from daily MODIS observations and then derive the median value of ρ for each grid from the 120 months of observation. As shown in Fig. 5e, at the subgrid level, the LWP and CDNC tend to be positively correlated almost over all tropical oceans. Mathematically, this is not surprising because as shown in Fig. 5b and c, the subgrid variability of r e is 1 order of magnitude smaller than that of LWP. Since CDNC is proportional to LWP 1 2 r −3 e according to Eq. (27), the subgrid variability of CDNC is mainly determined by the variability of LWP, leading to the positive correlation. Physically, the correlation can be explained by several mechanisms. For example, Wood et al. (2018) andO et al. (2018) found that a large amount of low-level water clouds over the stratocumulus-to-cumulus transition are "optically thin veil clouds". These clouds are usually associated with low LWP and low CDNC (therefore positive correlation) and probably caused by the strong precipitation scavenging process in the active cumulus. Note that our definition of ρ is the subgrid spatial correlation of LWP and CDNC. It may be different from the definition used in many aerosol indirect effect studies where the temporal correlation of monthly mean LWP and CDNC is more interesting.

Influence of subgrid variation of cloud water
As discussed in Sect. 2.2, most current GCMs only consider the impact of subgrid cloud water variation on autoconversion rate but ignore the impact of subgrid CDNC variation. To make our analysis relevant to the current GCMs, we first analyze E q in Eq. (20) based on observations. The impacts of subgrid CDNC variation (i.e., E N ) and its correlation with cloud water (i.e., E COV ) will be analyzed in the next section.
We derive E q using two approaches. First, we derive it from the observed LWP PDF based on Eq. (20). As such, we do not have to make any assumptions about the shape of LWP PDF, although solving the integration in Eq. (20) is time-consuming. In the second approach, we first derive the relative inverse relative variance v of LWP and then derive the enhancement factor by assuming the subgrid PDF to be either gamma or lognormal. This approach is more efficient, but it may be subject to error if the true PDF deviates from the assumed PDF shape. Figure 6a shows the annual mean enhancement factor E q in the tropical region derived based on Eq. (20) (i.e., the first approach) from 10 years of MODIS observations. Figure 6b and c show the annual mean enhancement factor E q derived by assuming the subgrid cloud water follows the lognormal (i.e., Eq. 14) and gamma distribution (i.e., Eq. 8), respectively. There are a couple of interesting and important points to note. First of all, similar to the gridmean quantities in Fig. 4, the enhancement factor E q also shows a clear Sc-to-Cu transition. Over the Sc decks, because clouds are more homogeneous ( v > 5), the enhancement factor E q is only around 1-2.5, while over the Cu regions, the more inhomogeneous clouds with v < 1 leads to a larger enhancement factor E q around 3-5. As mentioned previously, in the current CAM5.3, E q is assumed to be a constant of 3.2. While this value is within the observational range, it obviously cannot capture the Sc-to-Cu transition. In fact, the constant value 3.2 overestimates the E q over the Sc Figure 6. The annual mean factor for the KK2000 scheme due to subgrid variation of LWP computed (a) directly from observation, i.e., E q in Eq. (20), (b) from relative variance assuming lognormal PDF of LWP, i.e., E q in Eq. (14), and (c) from relative variance assuming the gamma PDF of LWP, i.e., E q in Eq. (8).
region and underestimates the E q over the Cu region, which could lead to unrealistic drizzle production in both regions and to consequential impacts on cloud water budget, radiation, and even aerosol indirect effects on the model. The second point to note is that the E q based on the lognormal PDF assumption in Fig. 6b agrees well with the results in Fig. 6a derived directly from the observation. In contrast, the E q based on the gamma PDF assumption in Fig. 6c tends to be smaller, especially in the Cu regions. This result seems to suggest that the lognormal distribution provides a better fit to the observed subgrid cloud water variation than the gamma distribution, which has rarely been noted and reported in the previous studies.
A flexible, cloud-regime-dependent E q could help improve the simulation of Sc-to-Cu transition in the GCM. If a GCM employs an advanced cloud parameterization scheme, such as CLUBB, that is able to provide regime-dependent information on subgrid cloud variation, i.e., v, then the enhancement factor E q could be diagnosed from v. However, most traditional cloud parameterization schemes do not provide information on subgrid cloud variation. In such cases, if one does not wish to use a constant E q but a varying regime-dependent scheme, then either v or E q need to be parameterized as a function of some grid-mean cloud prop-erties resolved by the GCM. In fact, several attempts have been made along this line. Based on the combination airborne in situ measurements and satellite remote sensing products, Boutle et al. (2014) parameterized the "fractional standard deviation" (which is equivalent to 1/ √ v in our definition) of liquid-phase cloud as a function of grid-mean cloud fraction. This scheme was later updated and tested in a host GCM in Hill et al. (2015) and was found to reduce the shortwave cloud radiative forcing biases in the model. In a recent study, Xie and Zhang (2015) derived the subgrid cloud variations from the ground-based observations from three DOE ARM sites, and then parameterized the inverse relative variance v as a function of the atmospheric stability. Figure 7a shows the variation of inverse relative variance v as a function of the grid-mean liquid-phase cloud fraction f liq . In general, the value of v increases with the increasing f liq , which is expected from the Sc-to-Cu increase of f liq in Fig. 4b and the Sc-to-Cu decrease of v in Fig. 5c. The v(f liq ) pattern in Fig. 7a is also consistent with the results reported in Wood and Hartmann (2006) and Lebsock et al. (2013). In the hope of obtaining a simple parameterization scheme for v(f liq ) that can be used in GCMs, we fit the median value of v as a simple third-order polynomial of f liq as follows: Figure 7. (a) The inverse relative variance v and (b) autoconversion enhancement factor due to LWP subgrid variability assuming lognormal PDF as a function of grid-mean liquid cloud fraction, where the solid line, dark shaded area, and light shaded area correspond to the median value, 25 %-75 % percentiles, and 10 %-90 % percentiles, respectively. The dotted lines correspond to simple third-order polynomial fitting. For reference, the v q (f liq ) parameterization scheme based on Boutle et al. (2014) is also plotted in panel (a).
For comparison, we also plotted the v(f liq ) parameterization developed by Boutle et al. (2014) for grid size of 100 km (dashed line in the figure). Apparently, comparing with Eq. (30), it generates slightly smaller v for a given the f liq . This difference is probably because our results are based on different data.
To test the performance of this simple parameterization, we first substitute the f liq from the MODIS daily mean level-3 product into the above equation and then use the resultant v to compute the enhancement factor E q . Unfortunately, the enhancement factor E q computed based on the parameterized v(f liq ) as shown in Fig. 8a substantially underestimates the observation-based results in Fig. 6, especially over the Cu regions. The deviation is probably because the relationship between E q and v is highly nonlinear (e.g., Eqs. 8 and 14), and therefore the above parameterization scheme that is designed for v is not able to capture the variability of E q . Based on this consideration, we tried an alternative approach. Instead of parameterization of v, we directly parameterize the enhancement factor E q as a function of f liq . Figure 7b shows the variation of E q as a function of f liq . As expected, E q generally decreases with increasing f liq . The median value of E q is fitted with the following third-order polynomial of f liq : As shown in Fig. 8b, the value of E q based on the above equation clearly agrees with the observation-based values in Fig. 6 better than that based on the parameterization of v(f liq ). The elimination of the middle step indeed improves the parameterization results. While this is encouraging, it should be kept in mind that Eq. (31) has very limited application; i.e., it is only useful for the autoconversion rate computation for a particular value of the autoconversion exponent beta, i.e., β q = 2.47. A good parameterization of v could be useful not only for autoconversion, but also for accretion and radiation computations. Another caution is that, if applied to a GCM, the performance of the E q (f liq ) parameterization in Eq. (31) will be dependent on the simulated accuracy of f liq in the model.

Influence of subgrid variance of CDNC
Now we will investigate the impacts of subgrid CDNC variation on the autoconversion rate simulation. For the moment, we will consider E N only. The impact of CDNC and cloud water correlation will be discussed in the next section. Similar to E q we first derive E N from the CDNC PDF based on Eq. (21). The annual mean result based on 10 years of MODIS observations is shown in Fig. 9a. There are several intriguing points to note. First of all, the value of E N is actually larger than E q in Fig. 9 such that we even have to use a different color scale for this plot. Secondly, the regions with escalated E N seem to coincide with the downwind regions of biomass burning aerosols (e.g., Gulf of Guinea, east coast of South Africa), air pollution (i.e., East China Sea), and, most interestingly, active volcanos (e.g., Kilauea, Hawaii, and Ambae, Vanuatu). We have also checked the seasonal variation of the E N and the results also support this observation. Another interesting feature to note is that, although the dust outflow regions, such as the tropical East Atlantic and Arabian Sea, have heavy aerosol loading, the value of E N there is only moderate. Figure 9b shows the value of E N computed based on Eq. (14) from the inverse relative variance of v, assuming that the subgrid CDNC follows a lognormal PDF. Although the overall pattern is consistent with Fig. 9a, the assumption of the lognormal PDF seems to underestimate E N . A closer examination indicates that the lognormal PDF tends to underestimate the population of clouds with a small CDNC and therefore underestimate the variance of CDNC as well as E N . We did not compute the E N based on the gamma distribution because of the singular value problem mentioned in Sect. 2.1. We could not find any previous observation-based study on the global pattern of the subgrid variation of CDNC and the corresponding E N . So, it is difficult for us to corroborate our results. On the one hand, the magnitude of E N is surprisingly large. As explained in Sect. 3, the CDNC is estimated based on Eq. (27) from the MODIS retrieval of COT and CER. Several previous studies have shown that the subpixellevel surface contamination, subpixel cloud inhomogeneity, and three-dimensional radiative transfer effects can cause significant errors in the MODIS CER retrievals, especially over broken cloud regions (Zhang and Platnick, 2011;Zhang et al., 2012Zhang et al., , 2016. Given the fact that the CDNC retrieval is highly sensitive to CER error as a result of N d ∼ r − 5 2 e , the influence of retrieval uncertainty on subgrid CDNC variation cannot be ruled out. On the other hand, the pattern of E N in Fig. 9a seems to suggest that there are some underlying physical mechanisms controlling the subgrid variation of CDNC, in which aerosols seem to play an important role. To achieve a better understanding, we analyzed the dependence of E N on liquid cloud fraction and grid-mean CDNC in Fig. 10, which reveals that E N has a stronger dependence on CDNC than cloud fraction. This result seems to indicate that the pattern of E N in Fig. 9 is largely determined by physical mechanisms rather than retrieval uncertainties. Interestingly, the largest E N is usually found when liquid cloud fraction is small and CDNC is large and decreases with decreasing CDNC and increasing cloud fraction. This pattern leads us to the following hypothesis: In the regions where aerosol is lim-ited, even weak updrafts can activate most cloud condensation nuclei (CCN). As a result, even if there is significant subgrid variation of turbulence at cloud base, the subgrid variation of CDNC remains small. In contrast, in regions where aerosol is abundant, the subgrid variation of turbulence becomes important. The subgrid variation of updraft leads to subgrid variation CDNC and thereby a large E N .
As far as we know, the results in Figs. 9 and 10 mark the first attempt based on satellite observations to unveil the global pattern of the subgrid variations of CDNC and investigate the consequential impacts on warm rain simulations in GCMs. Although obscured by satellite retrieval uncertainties, the results still provide valuable insights. First of all, the enhancement factor E N due to the subgrid variations of CDNC is non-negligible, even comparable to the effect of subgrid cloud water variation (i.e., E q ). Second, the global pattern of E N in Fig. 9 provides a valuable map for future studies.

The combined effect of subgrid variations of cloud water and CDNC
Finally, in this section we examine the combined effect of subgrid variations of cloud water and CDNC, as well as their correlation, on the autoconversion rate simulation. The annual mean combined enhancement factor E derived based on Eq. (17) from 10 years of MODIS COT and CER observations is shown in Fig. 11a. Comparing with the E q in Fig. 6 and E N in Fig. 9, the combined enhancement factor is generally larger. It is easy to see that in some regions (e.g., Gulf of Guinea, east coast of South Africa, and East China Sea), the combined enhancement factor E resembles E N , while in other regions (i.e., trade wind cumulus regions over open ocean) it resembles more of E q . Interestingly, because  both E q and E N are small over the Sc decks, those regions have the smallest combined enhancement factor E. As discussed in Sect. 2.2, only when the subgrid variation of cloud water is uncorrelated with the subgrid variation of CDNC can the combined enhancement factor E be decomposed into the simple product of E q and E N (i.e., Eq. 19). Figure 11b shows the annual mean value of the simple product E q · E N , without considering the correlation between cloud water and CDNC. Evidently, the simple product substantially overestimates the combined enhancement factor derived from the joint PDF of LWP and CDNC. This result can be explained by the mostly positive subgrid correlation between LWP and CDNC in Fig. 5e. As explained in Sect. 2.2, the positive cor-relation means that clouds with more water also tend to have more CDNC. The autoconversion rate of such configurations is lower than that when LWP and CDNC have no correlation. Together, the E q in Fig. 6, E N in Fig. 9, and the combined enhancement factor in Fig. 11 lead us to the following important conclusion. It is not sufficient to consider only the impact of subgrid variation of cloud water (i.e., E q ) on the autoconversion rate simulation. The influences of subgrid CDNC variation, as well as the correlation between cloud water and CDNC, must also be taken into account to avoid significant error.
Finally, the combined enhancement factor derived based on Eq. (23), assuming that the LWP and CDNC follow the bivariate lognormal distribution, is shown in Fig. 11c. Despite the tendency of overestimation, the result agrees reasonably well with that based on observed joint PDFs in Fig. 11a, clearly better than the simple product E q ·E N . This is encouraging as it suggests that the bivariate lognormal distribution can be used in the future to model the combined effect of cloud water and CDNC on autoconversion rate simulation in GCMs.

Summary and outlook
One of the difficulties in GCM simulation of the warm rain parameterization is how to account for the impact of subgrid variations of cloud properties, such as cloud water and CDCN, on nonlinear precipitation processes such as autoconversion. In practice, this impact is often treated by adding the enhancement factor term to the parameterization scheme. In this study, we derived the subgrid variations of liquidphase cloud properties over the tropical ocean using the satellite remote sensing products from MODIS and investigated the corresponding enhancement factors for parameterizations Figure 11. The combined enhancement factor derived (a) based on Eq. (17) from the observed joint PDF of LWP and CDNC, (b) assuming that subgrid variations of LWP and CDNC are uncorrelated, i.e., E q · E N only, and (c) based on Eq. (23) assuming that the subgrid LWP and CDNC following the bivariate lognormal distribution. of autoconversion rate. In comparison with previous work, our study is able to shed some new light on this problem in the following regards. 1. A theoretical framework is presented to explain the importance of the subgrid variation of CDNC and its correlation with cloud water on the autoconversion rate simulation in GCMs.
2. The wide spatial coverage of the level-3 MODIS product enables us to depict a detailed quantitative picture of the enhancement factor E q , which shows a clear cloud regime dependence, i.e., a Sc-to-Cu increase. The constant E q = 3.2 used in the current CAM5.3 model overestimates and underestimates the observed E q in the Sc and Cu regions, respectively.
3. The E q based on the lognormal PDF assumption performs significantly better than that based on the gamma PDF assumption. A simple parameterization scheme is provided to relate E q to the grid-mean liquid cloud fraction, which can be readily used in GCMs.
4. For the first time, the enhancement factor E N due to the subgrid variation of CDNC is derived from satellite observation, and the results reveal several regions downwind of biomass burning aerosols (e.g., Gulf of Guinea, east coast of South Africa), air pollution (i.e., East China Sea), and active volcanos (e.g., Kilauea, Hawaii, and Ambae, Vanuatu). The largest E N is usually found where CDNC is large and liquid cloud fraction is small and decreases with decreasing CDNC and increasing cloud fraction.
5. MODIS observations suggest that the subgrid LWP and CDNC are mostly positively correlated. As a result, the combined enhancement factor is significantly smaller than the simple product of E q · E N (i.e., assuming no correlation). The combined enhancement factor derived assuming LWP and CDNC follow the bivariate lognormal distribution agrees with the observation-based results reasonably well.
As noted in the previous sections, this study has several important limitations, most of which are a result of using the level-3 MODIS observations. The fixed 1 • × 1 • spatial resolution of the MODIS level-3 product makes it impossible for us to investigate the scale dependence of subgrid cloud variation. Similar to previous studies, we have to make several assumptions when estimating the CDNC from the level-3 MODIS product. Furthermore, the retrieval uncertainties associated with the optically thin clouds in the MODIS product pose a challenging obstacle for the quantification of subgrid cloud property variations and the corresponding enhancement factors. These limitations have to be addressed using additional independent observations from, for example, ground-based remote sensing products and/or in situ measurements from airborne field campaigns. Recently, a few novel methods have been developed to provide certain information on the subgrid cloud property variations to the host GCM. Most noticeable examples are the superparameterization method (a.k.a. multiscale modeling framework)  and the PDF-based higher-order turbulence closure methods, e.g., Cloud Layer Unified By Binormals (CLUBB) (Golaz et al., 2002a;Guo et al., 2015;Larson et al., 2002) and Eddy-Diffusivity Mass Flux (EDMF) (Sušelj et al., 2013). The subgrid cloud property variations derived in this study provide the valuable observational basis for the evaluation and improvement of these schemes. Data availability. The MODIS level-3 daily averaged products used in this study are available from the Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) of NASA. All results from the analysis can be made available upon request to the corresponding author.