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 (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.

1 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

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial q_{\mathrm{r}}}{\partial t}}=C{\left(q_{\mathrm{c}}\right)}^{{\mathit{\beta }}_q}{\left(N_{\mathrm{c}}\right)}^{{\mathit{\beta }}_N},\end{array}$$

where $\frac{\partial q_{\mathrm{r}}}{\partial t}$ is the rain water tendency due to the autoconversion process, q_c has the unit of kg kg⁻¹, and N_c the unit of cm⁻³. The three parameters C=1350, β_q=2.47, and ${\mathit{\beta }}_N=-\mathrm{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 process-level 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 grid-mean 〈x〉, i.e., 〈f(x)〉≠f(〈x〉). Mathematically, if f(x) is convex, then f(〈x〉)<〈f(x)〉 (Larson and Griffin, 2013; 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 $\langle f\left(x\right)\rangle =E\cdot f(\langle x\rangle )$. 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 (Barker et al., 1996; 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 large-eddy 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 near-zero cloud water. In addition to CLUBB, the so-called super-parameterization (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 liquid-phase 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.

2 Theoretical background

2.1 Theoretical distributions to describe subgrid cloud property variations

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):

$$\begin{array}{}\text{(2)}& {\displaystyle }{P}_{\mathrm{G}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left(v\right)}}{\mathit{\alpha }}^{\mathit{\nu }}{x}^{\mathit{\nu }-\mathrm{1}}\exp \left(-\mathit{\alpha }x\right),\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\displaystyle }\langle x\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}x{P}_{\mathrm{G}}\left(x\right)\mathrm{d}x={\displaystyle \frac{v}{\mathit{\alpha }}}\end{array}$$

and variance given by

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{Var}\left(x\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(x-\langle x\rangle \right)}^{\mathrm{2}}{P}_{\mathrm{G}}\left(x\right)\mathrm{d}x={\displaystyle \frac{v}{{\mathit{\alpha }}^{\mathrm{2}}}}.\end{array}$$

It follows from Eqs. (3) and (4) that the so-called inverse relative variance is

$$\begin{array}{}\text{(5)}& {\displaystyle }v={\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}={\displaystyle \frac{\langle x{\rangle }^{\mathrm{2}}}{\mathrm{Var}\left(x\right)}},\end{array}$$

where $\mathit{\eta }=\frac{\mathrm{Var}\left(x\right)}{\langle x{\rangle }^{\mathrm{2}}}$ is the relative variance. If x follows the gamma distribution for a physical process M(x) that is a power function of x,

$$\begin{array}{}\text{(6)}& {\displaystyle }M\left(x\right)=K{x}^{\mathit{\beta }},\end{array}$$

then the expected value 〈M(x)〉 is given by

$$\begin{array}{ll}{\displaystyle }\langle M\left(x\right){\rangle }_G& {\displaystyle }=K\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^{\mathit{\beta }}{P}_{\mathrm{G}}\left(x\right)\mathrm{d}x\\ \text{(7)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{\Gamma }(v+\mathit{\beta })}{\mathrm{\Gamma }\left(v\right)v^{\mathit{\beta }}}}K\langle \left(x\right){\rangle }^{\mathit{\beta }},\mathit{\beta }>-v.\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }E\left(P_{\mathrm{G}},v,\mathit{\beta }\right)& {\displaystyle }={\displaystyle \frac{\langle K{x}^{\mathit{\beta }}\rangle }{K\langle x{\rangle }^{\mathit{\beta }}}}={\displaystyle \frac{\mathrm{1}}{\langle x{\rangle }^{\mathit{\beta }}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^{\mathit{\beta }}{P}_{\mathrm{G}}\left(x\right)\mathrm{d}x\\ \text{(8)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{\Gamma }(v+\mathit{\beta })}{\mathrm{\Gamma }\left(v\right)v^{\mathit{\beta }}}}.\end{array}$$

The PDF of a lognormal distribution is given as follows (Larson and Griffin, 2013; Lebsock et al., 2013):

$$\begin{array}{}\text{(9)}& {\displaystyle }{P}_{\mathrm{L}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}x\mathit{\sigma }}}\exp \left(-{\displaystyle \frac{{\left(\ln x-\mathit{\mu }\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where μ=〈lnx〉 and σ²=Var(lnx) correspond to the mean and variance of lnx, respectively. The mean value of the lognormal distribution is given by

$$\begin{array}{}\text{(10)}& {\displaystyle }\langle x\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}x{P}_{\mathrm{L}}\left(x\right)\mathrm{d}x=e^{\mathit{\mu }+\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}}\end{array}$$

and the variance by

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathrm{Var}\left(x\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(x-\langle x\rangle \right)}^{\mathrm{2}}{P}_{\mathrm{L}}\left(x\right)\mathrm{d}x=e^{\mathrm{2}\mathit{\mu }+{\mathit{\sigma }}^{\mathrm{2}}}\left(e^{{\mathit{\sigma }}^{\mathrm{2}}}-\mathrm{1}\right).\end{array}$$

It follows from Eqs. (10) and (11) that the inverse relative variance can be derived from the following equation:

$$\begin{array}{}\text{(12)}& {\displaystyle }{e}^{{\mathit{\sigma }}^{\mathrm{2}}}=\mathrm{1}+{\displaystyle \frac{\mathrm{Var}\left(x\right)}{\langle x{\rangle }^{\mathrm{2}}}}=\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{v}}.\end{array}$$

If x follows the lognormal distribution, the expected value of 〈M(x)〉 is

$$\begin{array}{}\text{(13)}& {\displaystyle }\langle M\left(x\right){\rangle }_L=K\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^{\mathit{\beta }}{P}_{\mathrm{L}}\left(x\right)\mathrm{d}x={\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{v}}\right)}^{\frac{{\mathit{\beta }}^{\mathrm{2}}-\mathit{\beta }}{\mathrm{2}}}K\langle x{\rangle }^{\mathit{\beta }}.\end{array}$$

Evidently, the corresponding enhancement factor is given by

$$\begin{array}{}\text{(14)}& {\displaystyle }E\left(P_{\mathrm{L}},v,\mathit{\beta }\right)={\displaystyle \frac{\langle K{x}^{\mathit{\beta }}\rangle }{K\langle x{\rangle }^{\mathit{\beta }}}}={\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{v}}\right)}^{\frac{{\mathit{\beta }}^{\mathrm{2}}-\mathit{\beta }}{\mathrm{2}}}.\end{array}$$

Note that Eqs. (7) and (8) are only valid when $\mathit{\beta }>-v$ because gamma function Γ(v+β) can run into singular values when $v+\mathit{\beta }<\mathrm{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.

Figure 1 (a) The probability density function (PDF) and cumulative distribution function (CDF) of cloud LWC (q_c) that follow the gamma (dashed) and lognormal (solid) distribution. For the both distributions, 〈q_c〉=0.5g kg⁻¹ and v_q=3.0. (b) The PDF and CDF of autoconversion rate computed based on the KK2000 scheme in Eq. (15) and the PDF of q_c. In the computation, the N_c is kept at a constant of 50 cm⁻¹. (c) The PDF and CDF of N_c that follow the gamma (dashed) and lognormal (solid) distribution. For the both distributions, 〈N_c〉=50 cm⁻³ and v_N=5.0. (d) The PDF and CDF of the autoconversion rate computed based on the KK2000 scheme in Eq. (15) and the PDF of N_c. The q_c is kept at 0.5 g kg⁻¹ in the computation.

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⁻¹ 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⁻¹, 〈N_c〉=50 cm⁻³, 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⁻¹) 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_{\mathrm{c}}^{\mathrm{2.47}}$.

Figure 2 Enhancement factors based on lognormal E(P_L, β) and gamma E(P_G, β) subgrid PDFs for different β as a function of the inverse relative variance v.

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.

2.2 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:

$$\begin{array}{}\text{(15)}& {\displaystyle }〈{\displaystyle \frac{\partial q_{\mathrm{r}}}{\partial t}}〉=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}C{\left(q_{\mathrm{c}}\right)}^{{\mathit{\beta }}_q}{\left(N_{\mathrm{c}}\right)}^{{\mathit{\beta }}_N}P\left(q_{\mathrm{c}},N_{\mathrm{c}}\right)\mathrm{d}{q}_{\mathrm{c}}\mathrm{d}{N}_{\mathrm{c}},\end{array}$$

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 following equation:

$$\begin{array}{}\text{(16)}& {\displaystyle }〈{\displaystyle \frac{\partial q_{\mathrm{r}}}{\partial t}}〉=E\cdot C{\left(\langle q_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_q}{\left(\langle N_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_N},\end{array}$$

where E is the enhancement factor defined as

$$\begin{array}{}\text{(17)}& {\displaystyle }E={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\left(q_{\mathrm{c}}\right)}^{{\mathit{\beta }}_q}{\left(N_{\mathrm{c}}\right)}^{{\mathit{\beta }}_N}P\left(q_{\mathrm{c}},N_{\mathrm{c}}\right)\mathrm{d}{q}_{\mathrm{c}}\mathrm{d}{N}_{\mathrm{c}}}{{\left(\langle q_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_q}{\left(\langle N_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_N}}}.\end{array}$$

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\left(q_{\mathrm{c}},N_{\mathrm{c}}\right)=P\left(q_{\mathrm{c}}\right)P\left(N_{\mathrm{c}}\right)$, where P(q_c) and P(N_c) are the PDF of the subgrid q_c and N_c. Consequently, Eq. (15) reduces to

$$\begin{array}{}\text{(18)}& {\displaystyle }〈{\displaystyle \frac{\partial q_{\mathrm{r}}}{\partial t}}〉=C\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(q_{\mathrm{c}}\right)}^{{\mathit{\beta }}_q}P\left(q_{\mathrm{c}}\right)\mathrm{d}{q}_{\mathrm{c}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(N_{\mathrm{c}}\right)}^{{\mathit{\beta }}_N}P\left(N_{\mathrm{c}}\right)\mathrm{d}{N}_{\mathrm{c}}\end{array}$$

and Eq. (17) to

$$\begin{array}{}\text{(19)}& {\displaystyle }E=E_q\cdot E_N,\end{array}$$

where E_q is the enhancement factor due to the subgrid variation of cloud water, which has the form

$$\begin{array}{}\text{(20)}& {\displaystyle }{E}_q={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\left(q_{\mathrm{c}}\right)}^{{\mathit{\beta }}_q}P\left(q_{\mathrm{c}}\right)\mathrm{d}{q}_{\mathrm{c}}}{{\left(\langle q_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_q}}},\end{array}$$

and the E_N is the enhancement factor due to the subgrid variation of CDNC, which has the form

$$\begin{array}{}\text{(21)}& {\displaystyle }{E}_N={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\left(N_{\mathrm{c}}\right)}^{{\mathit{\beta }}_N}P\left(N_{\mathrm{c}}\right)\mathrm{d}{N}_{\mathrm{c}}}{{\left(\langle N_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_N}}}.\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }P\left(q_{\mathrm{c}},N_{\mathrm{c}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{q}_{\mathrm{c}}{N}_{\mathrm{c}}{\mathit{\sigma }}_q{\mathit{\sigma }}_N\sqrt{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}}\exp \left(-{\displaystyle \frac{\mathit{\zeta }}{\mathrm{2}}}\right)\\ {\displaystyle }& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\rho }}^{\mathrm{2}}}}\left[{\left({\displaystyle \frac{\ln q_{\mathrm{c}}-{\mathit{\mu }}_q}{{\mathit{\sigma }}_q}}\right)}^{\mathrm{2}}-\mathrm{2}\mathit{\rho }\left({\displaystyle \frac{\ln q_{\mathrm{c}}-{\mathit{\mu }}_q}{{\mathit{\sigma }}_q}}\right)\left({\displaystyle \frac{\ln N_{\mathrm{c}}-{\mathit{\mu }}_N}{{\mathit{\sigma }}_N}}\right)\right.\\ \text{(22)}& {\displaystyle }& {\displaystyle }\left.+{\left({\displaystyle \frac{\ln N_{\mathrm{c}}-{\mathit{\mu }}_N}{{\mathit{\sigma }}_N}}\right)}^{\mathrm{2}}\right],\end{array}$$

where ρ is the correlation coefficient between q_c and N_c (Larson and Griffin, 2013; 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:

$$\begin{array}{ll}{\displaystyle }E& {\displaystyle }=E_q\left(P_{\mathrm{L}},v_q,{\mathit{\beta }}_q\right)\cdot E_N\left(P_{\mathrm{L}},v_N,{\mathit{\beta }}_N\right)\\ \text{(23)}& {\displaystyle }& {\displaystyle }\cdot E_{\mathrm{COV}}\left(\mathit{\rho },{\mathit{\beta }}_q,{\mathit{\beta }}_N{v}_q,v_N\right),\end{array}$$

where $E_q\left(P_{\mathrm{L}},v_q,{\mathit{\beta }}_q\right)={\left(\mathrm{1}+\frac{\mathrm{1}}{v_q}\right)}^{\frac{{\mathit{\beta }}_q^{\mathrm{2}}-{\mathit{\beta }}_q}{\mathrm{2}}}$ and E_N (P_L, v_N, β_N) = ${\left(\mathrm{1}+\frac{\mathrm{1}}{v_N}\right)}^{\frac{{\mathit{\beta }}_N^{\mathrm{2}}-{\mathit{\beta }}_N}{\mathrm{2}}}$ correspond to the impacts of subgrid q_c and N_c variance, respectively (i.e., Eq. 14), and

$$\begin{array}{}\text{(24)}& {\displaystyle }{E}_{\mathrm{COV}}\left(\mathit{\rho },{\mathit{\beta }}_q,{\mathit{\beta }}_N,v_q,v_N\right)=\exp \left(\mathit{\rho }{\mathit{\beta }}_q{\mathit{\beta }}_N{\mathit{\sigma }}_q{\mathit{\sigma }}_N\right)\end{array}$$

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., $\mathit{\rho }=\mathrm{0},E_{\mathrm{COV}}=\mathrm{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., 2015, 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, 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., ${\int }_{\mathrm{0}}^{\mathrm{\infty }}{\left(N_{\mathrm{c}}\right)}^{{\mathit{\beta }}_N}P\left(N_{\mathrm{c}}\right)\mathrm{d}{N}_{\mathrm{c}})$ would be different from that computed based on the mean of N_c (i.e., ${\left(\langle N_{\mathrm{c}}\rangle \right)}^{{\mathit{\beta }}_N})$. The autoconversion enhancement factor based on the lognormal PDF E(P_L,β) for ${\mathit{\beta }}_N=-\mathrm{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 ${\mathit{\beta }}_N=-\mathrm{1.79}$ is actually larger than that for β_q=2.47 because of the formula of the exponent in Eq. (14) (i.e., $\frac{{\mathit{\beta }}^{\mathrm{2}}-\mathit{\beta }}{\mathrm{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\cdot E_{\mathrm{COV}}>\mathrm{1}$. After some manipulation, we can show that if β_N<0 and σ_N>0, then

$$\begin{array}{}\text{(25)}& {\displaystyle }{E}_N\cdot E_{\mathrm{COV}}>\mathrm{1},\mathrm{if}\mathit{\rho }<{\displaystyle \frac{{\mathit{\sigma }}_N}{{\mathit{\sigma }}_q}}\cdot {\displaystyle \frac{\left(\mathrm{1}-{\mathit{\beta }}_N\right)}{\mathrm{2}{\mathit{\beta }}_q}}.\end{array}$$

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.

3 Data and methodology

To derive the above-mentioned enhancement factors, we will use 10 years (2007–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 pixel-level (i.e., level-2) retrievals. As summarized in Platnick et al. (2003, 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 grid-mean values (〈x〉) and subgrid standard deviations (σ_x). The inverse relative variance v can then be derived from Eq. (5), i.e., $v=\langle x{\rangle }^{\mathrm{2}}/{\mathit{\sigma }}_x^{\mathrm{2}}$. 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., 2012, 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.

Figure 3 The (a) LWP and (b) CDNC as a function of COT and CER. (c) An example of the COT–CER joint histogram observed by Aqua MODIS on 9 January 2007 at 1^∘ S and 1^∘ W.

Given the COT and CER retrieval, the operational MODIS product estimates the LWP of clouds using

$$\begin{array}{}\text{(26)}& {\displaystyle }\mathrm{LWP}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{w}}\mathrm{COT}\cdot \mathrm{CER},\end{array}$$

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=\langle x{\rangle }^{\mathrm{2}}/{\mathit{\sigma }}_x^{\mathrm{2}}$. 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 CloudSat 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

$$\begin{array}{ll}{\displaystyle }{N}_{\mathrm{c}}\left(\mathit{\tau },r_{\mathrm{e}}\right)& {\displaystyle }={\displaystyle \frac{\sqrt{\mathrm{5}}}{\mathrm{2}\mathit{\pi }k}}{\displaystyle \frac{\sqrt{f_{\mathrm{ad}}{\mathrm{\Gamma }}_{\mathrm{w}}}}{\sqrt{{\mathit{\rho }}_{\mathrm{w}}{Q}_{\mathrm{e}}}}}{\mathit{\tau }}^{\frac{\mathrm{1}}{\mathrm{2}}}{r}_{\mathrm{e}}^{-\frac{\mathrm{5}}{\mathrm{2}}}\\ \text{(27)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\sqrt{\mathrm{15}}}{\mathrm{2}\mathit{\pi }k}}{\displaystyle \frac{\sqrt{f_{\mathrm{ad}}{\mathrm{\Gamma }}_{\mathrm{w}}}}{{\mathit{\rho }}_{\mathrm{w}}{\sqrt{\mathrm{2}Q}}_{\mathrm{e}}}}{\mathrm{LWP}}^{\frac{\mathrm{1}}{\mathrm{2}}}{r}_{\mathrm{e}}^{-\mathrm{3}},\end{array}$$

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 volume-equivalent 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 time-consuming, 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., 2017, 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:

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }\langle x\rangle =\int \int x\left(\mathit{\tau },r_{\mathrm{e}}\right)P\left(\mathit{\tau },r_{\mathrm{e}}\right)\mathrm{d}\mathit{\tau }\mathrm{d}{r}_{\mathrm{e}},\text{(29)}& {\displaystyle }& {\displaystyle }\mathrm{Var}\left(x\right)=\int \int {\left(x\left(\mathit{\tau },r_{\mathrm{e}}\right)-\langle N_{\mathrm{c}}\rangle \right)}^{\mathrm{2}}P\left(\mathit{\tau },r_{\mathrm{e}}\right)\mathrm{d}\mathit{\tau }\mathrm{d}{r}_{\mathrm{e}},\end{array}$$

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.

Figure 4 Ten-year (2007–2016) averaged annual mean (a) total cloud fraction, (b) liquid cloud fraction, (c) cloud optical thickness, (d) cloud effective radius retrieved from the 3.7 µm band, (e) cloud water path, and (f) cloud droplet concentration retrievals from Aqua MODIS over the tropical (30^∘ S–30^∘ N) oceans. All quantities are in-cloud means that are averaged over the cloudy part of the grid only.

4 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), liquid-phase 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 prominent 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).

Figure 5 Median value of the inverse relative variance (i.e., $v=\langle x{\rangle }^{\mathrm{2}}/\mathrm{Var}\left(x\right)$) for (a) COT, (b) CER, (c) LWP, and (d) CDNC and (e) median value of the correlation coefficient between LWP and CDNC derived from 10 years of MODIS observations. Note that the color scale of CER is different from others'.

Switching the focus now from grid-mean values to subgrid variability, we will show the grid-level inverse relative variances $v=\langle x{\rangle }^{\mathrm{2}}/\mathrm{Var}\left(x\right))$ 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 $\stackrel{\mathrm{̃}}{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 $\stackrel{\mathrm{̃}}{v}$ values of all four sets of cloud properties (i.e., COT, CER, LWP, and CDNC) all exhibit a clear and similar Sc-to-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 $\stackrel{\mathrm{̃}}{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 $\stackrel{\mathrm{̃}}{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_{\mathrm{c}}\sim r_{\mathrm{e}}^{-\frac{\mathrm{5}}{\mathrm{2}}}$) 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 $\stackrel{\mathrm{̃}}{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 $\stackrel{\mathrm{̃}}{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 $\mathit{\chi }=\exp \left(\ln \langle \mathit{\tau }\rangle \right)/\langle \mathit{\tau }\rangle$), 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 $\stackrel{\mathrm{̃}}{\mathit{\rho }}$. Similar to the derivation of median $\stackrel{\mathrm{̃}}{v}$, we first compute the monthly mean ρ from daily MODIS observations and then derive the median value of $\stackrel{\mathrm{̃}}{\mathit{\rho }}$ 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 ${\mathrm{LWP}}^{\frac{\mathrm{1}}{\mathrm{2}}}{r}_{\mathrm{e}}^{-\mathrm{3}}$ 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) and O 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.

5 Implications for warm rain simulations in GCM

5.1 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.

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).

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 grid-mean 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 ($\stackrel{\mathrm{̃}}{v}>\mathrm{5}$), the enhancement factor E_q is only around 1–2.5, while over the Cu regions, the more inhomogeneous clouds with $\stackrel{\mathrm{̃}}{v}<\mathrm{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 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 properties 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 $\mathrm{1}/\sqrt{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 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).

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:

$$\begin{array}{ll}{\displaystyle }v\left(f_{\mathrm{liq}}\right)& {\displaystyle }=\mathrm{2.38}-\mathrm{4.95}{f}_{\mathrm{liq}}+\mathrm{8.74}{f}_{\mathrm{liq}}^{\mathrm{2}}\\ \text{(30)}& {\displaystyle }& {\displaystyle }-\mathrm{0.49}{f}_{\mathrm{liq}}^{\mathrm{3}},f_{\mathrm{liq}}\in [\mathrm{0},\mathrm{1}].\end{array}$$

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.

Figure 8 Annual mean value of the enhancement factor E_N computed based on the (a) $v\left(f_{\mathrm{liq}}\right)=\mathrm{2.38}-\mathrm{4.95}{f}_{\mathrm{liq}}+\mathrm{8.74}{f}_{\mathrm{liq}}^{\mathrm{2}}-\mathrm{0.49}{f}_{\mathrm{liq}}^{\mathrm{3}}$ parameterization scheme in Eq. (30) and (b) $E_q\left(f_{\mathrm{liq}}\right)=\mathrm{2.72}+\mathrm{7.33}{f}_{\mathrm{liq}}-\mathrm{19.17}{f}_{\mathrm{liq}}^{\mathrm{2}}+\mathrm{10.69}{f}_{\mathrm{liq}}^{\mathrm{3}}$ parameterization scheme in Eq. (31).

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:

$$\begin{array}{ll}{\displaystyle }{E}_q\left(f_{\mathrm{liq}}\right)& {\displaystyle }=\mathrm{2.72}+\mathrm{7.33}{f}_{\mathrm{liq}}-\mathrm{19.17}{f}_{\mathrm{liq}}^{\mathrm{2}}\\ \text{(31)}& {\displaystyle }& {\displaystyle }+\mathrm{10.69}{f}_{\mathrm{liq}}^{\mathrm{3}},f_{\mathrm{liq}}\in [\mathrm{0},\mathrm{1}].\end{array}$$

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.

5.2 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.

Figure 9 Annual mean value of the enhancement factor E_N derived from (a) observation based on Eq. (21) and (b) from Eq. (14) assuming lognormal subgrid CDNC distribution.

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 subpixel-level 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., 2012, 2016). Given the fact that the CDNC retrieval is highly sensitive to CER error as a result of $N_{\mathrm{d}}\sim r_{\mathrm{e}}^{-\frac{\mathrm{5}}{\mathrm{2}}}$, 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 limited, 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.

Figure 10 Dependence of E_N on f_liq and N_d. The color map corresponds to the mean value of E_N for a given N_d and f_liq bin. The white contour lines correspond to the relative sampling frequency of N_d and f_liq bins (i.e., the most frequently observed combination is N_d∼50 cm⁻³ and f_liq∼0.1).

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.

5.3 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 correlation 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.

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.

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.

6 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 liquid-phase cloud properties over the tropical ocean using the satellite remote sensing products from MODIS and investigated the corresponding enhancement factors for parameterizations 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 super-parameterization method (a.k.a. multiscale modeling framework) (Wang et al., 2015) 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.