Improving the Treatment of Subgrid Cloud Variability in Warm Rain Simulation in CESM2

Representing subgrid variability of cloud properties has always been a challenge in global climate models (GCMs). In many cloud microphysics schemes, the warm rain non‐linear process rates calculated based on grid‐mean cloud properties are usually scaled by an enhancement factor (EF) to account for the effects of subgrid cloud variability. In our study, we find that the EF derived from Cloud Layers Unified by Binormals in Community Atmosphere Model version 6 (CAM6) is severely overestimated in most of the cloudy oceanic areas, which leads to strong overestimation of the autoconversion rate. We improve the EF in warm rain simulation by developing a new formula for in‐cloud subgrid cloud water variance. With the updated subgrid cloud water variance and EF treatment, the liquid cloud fraction (LCF) and cloud optical thickness (COT) increases noticeably for marine stratocumulus, and the shortwave cloud forcing (SWCF) matches better with observations. The updated formula improves the relationship between autoconversion rate and cloud droplet number concentration, and it decreases the sensitivity of autoconversion rate to aerosols. The sensitivity of liquid water path to aerosols decreases noticeably and is in better agreement with that in MODIS. Although the sensitivity of COT is similar to that in MODIS, CAM6 underestimates the sensitivity of grid‐mean SWCF to aerosols because of the underestimation in the sensitivities of LCF and in‐cloud SWCF. Our results indicate the importance of representing reasonable subgrid cloud variability in the simulation of cloud properties and aerosol‐cloud interaction in GCMs.

The accuracy of subgrid parameterization directly affects the calculation of the autoconversion rate. Autoconversion rate ( ) auto is usually parameterized as a function of in-cloud LWC and CDNC. Khairoutdinov and Kogan provided a particular formula of autoconversion rate based on large-eddy simulation with bin microphysics that can simulate the process-level physics (Khairoutdinov & Kogan, 2000: where α = 2.47, β = −1.79, the units of ( ) auto , LWC and CDNC are kgkg −1 s −1 , kgkg −1 , and cm −3 . This formula is widely used in GCMs. However, only mean variables, LWC and CDNC , are available in GCMs, and so they have to be used to calculate the mean autoconversion rate. If a function ( ) is nonlinearly dependent on a variable , the untuned ( ) based on the mean value is not equal to the desired mean value ( ) that could have been computed if the subgrid values of were known (Pincus & Klein, 2000). In the parameterization of autoconversion, this inequality is often accounted for through tuning the autoconversion rate ( ) by an enhancement factor (EF) that is, ( ) = ⋅ ( ) (Gettelman et al., 2008a(Gettelman et al., , 2008b. EF is determined primarily by two factors: the nonlinearity of its function and the subgrid probability density function (PDF) of LWC. EF is larger when the function is more nonlinear or when the LWC is more inhomogeneous in the cloudy part of the grid point.
Because it is challenging for GCMs to simulate the subgrid variability of cloud properties, the previous versions of Community Atmosphere Model (CAM) usually do not provide the information of subgrid cloud water variance (Meehl et al., 2013). As a result, a simple constant EF = 3.2 is assumed, which is based on an observational study by Barker et al. (1996) that found the typical value of EF around 3.2 for marine stratocumulus clouds in satellite observations. Obviously, a constant EF = 3.2 cannot capture the variation of EF in different cloud regimes (Boutle et al., 2014;Zhang et al., 2019). In the up-to-date version CAM6, which is the atmospheric component of the Community Earth System Model version 2 (CESM2; Hurrell et al., 2013), the advanced macrophysics scheme Cloud Layers Unified by Binormals (CLUBB) has been implemented, which can provide the subgrid variances of cloud water diagnosed from the joint PDFs. The implementation of CLUBB makes it possible to dynamically diagnose EF through subgrid cloud water variance.
However, challenges remain in CAM6 in the simulations of cloud properties and their subgrid variance. Although CLUBB improves the simulation of clouds especially in the stratocumulus to cumulus (Sc-to-Cu) transition areas, it still underestimates in-cloud water and radiative forcing in MBL clouds (Bogenschutz et al., 2013;Z. Guo et al., 2014;Kubar et al., 2015;Zheng et al., 2016). It also has been revealed that a preliminary version of CLUBB implemented in CAM5 overestimates the EF in trade wind cumulus cloud regions, resulting in the overlarge autoconversion rate and overly fast consumption of cloud water (Song et al., 2018), which may partly contribute to the aforementioned problems. It is still unclear what might lead to the overestimation in EF. This indicates further work is needed to better understand how EF is treated and to improve EF treatment in climate models.
The EF treatment is also important for simulating cloud microphysical processes and ACI. Cloud properties typically have larger sensitivity to aerosols in GCMs than satellite observations (Lohmann & Lesins, 2002;Quaas & Boucher, 2005;Quaas et al., 2008Quaas et al., , 2010Sato et al., 2018), which is highly related to EF and autoconversion rate. It is important to examine how EF treatments affect ACI simulations in climate models.
In this study, we improve the treatment of subgrid cloud variance and EF in CAM6-CLUBB through a new formula. We will focus on the improvements in the liquid-phase clouds over the ocean. In Section 2, the calculation WANG ET AL.

10.1029/2022MS003103
3 of 20 of subgrid cloud water variance and EF, the satellite data, and the model settings are introduced. In Section 3, we use satellite observations to evaluate model performance in simulating subgrid cloud water variance, EF, cloud properties, and ACI over oceanic areas. The summary and discussion are included in Section 4.

Model Description
CESM2 is a state-of-the-art climate model consisting of atmosphere, land, ocean, and sea-ice components that exchange information and fluxes from each component via a coupler. The atmospheric component of CESM2 is CAM6. The primary change of CAM6 compared to previous versions is the inclusion of the unified turbulence scheme, CLUBB, which is a partial third-order turbulence closure scheme that utilizes a double-Gaussian joint PDF to model the variance and covariance of subgrid vertical velocity, temperature, and moisture Larson et al., 2002). It directly replaces the separate shallow-convection, boundary layer, and grid-scale condensation schemes present in CAM5 (Bogenschutz et al., 2013).
The prognostic two-moment stratiform cloud microphysics scheme by Morrison and Gettelman (Gettelman et al., 2015, hereafter MG2) is used in CAM6 for both the stratiform and shallow convection clouds. The main unique aspect of MG2 relative to the other schemes developed for GCMs with prognostic precipitation is the prognostic two-moment approach for all hydrometeor categories . MG2 uses the grid-mean value of LWC and CDNC together with the variance of subgrid cloud water diagnosed from CLUBB to calculate the microphysical processes rates.
It is important to note that the KK2000 formula of autoconversion rate in MG2 has been adjusted as follows in the default settings of CAM6: The units of the variables are same as that in Equation 1. EF is the enhancement factor of autoconversion rate. The exponent of CDNC is reduced (−1.1) following data analysis with VAMOS Ocean-Cloud-Atmosphere-Land Study data, and the prefactor is tuned (0.01) to get the cloud liquid water amount similar to previous model versions (Gettelman et al., 2019). Also, there is restriction on autoconversion rates in MG2. If cloud water consumption by the precipitation is larger than the cloud water amount within one microphysics time step (5 min), MG2 will adjust the processes rates including autoconversion rate to the value which removes all of the cloud water (Gettelman et al., 2008a(Gettelman et al., , 2008b.

The Calculation of Subgrid Cloud Water Variance and EF
MG2 assumes that the PDF of in-cloud cloud water ( (LWC)) follows a gamma distribution function (Barker, 1996;Oreopoulos & Davies, 1998): ν is the inverse relative variance of cloud water calculated by CLUBB; LWC is the mean in-cloud cloud water content in a grid point; Var(LWC) is the subgrid variance of in-cloud cloud water; = LWC ; Γ is the gamma function. Then the EF of KK2000 in MG2 is (Gettelman et al., 2008a(Gettelman et al., , 2008bZhang et al., 2019): Increased ν leads to decreased EF, and vice-versa.
In the default settings of CAM6, CLUBB uses the grid-mean cloud water content (LWC grid ) and double-Gaussian joint PDFs including both clear-sky and cloudy part to calculate the grid-mean subgrid cloud water variance (Var ( LWCgrid ) ): 4 of 20 χ is an "extended" liquid water mixing ratio, which represents liquid water mixing ratio when χ ≥ 0; H is the Heaviside step function; ( ) is the ith PDF of (two in total for the double-Gaussian PDF scheme in CLUBB) and is applied to both clear-sky and cloudy part; is the weight of the ith PDF, and 1 + 2 = 1 ; = 1 1 + 2 2 is the cloud fraction in a grid point; 1 and 2 are the cloud fraction derived by each CLUBB PDF. Then the inverse relative variance of cloud water is calculated as: As is shown later in Section 3, this calculation strongly overestimates the subgrid cloud water variance especially over the Cu regions where cloud fraction is small, because grid-mean values are used in calculating ν. For example, for a double delta function PDF, that is, a uniform cloud and a uniform environment, the within-cloud variance is zero, by definition. However, the grid-mean variance might be very large if the in-cloud liquid is large according to Equation 6.
Here we update the formula of ν by using the in-cloud values instead. We first derive in-cloud cloud water variance through the combination of two CLUBB Gaussian PDF components: The specific derivations of Equation 6 and Equation 8 are in Appendix A.
H. Guo et al. (2014) gives another formula of the in-cloud cloud water variance: They obtain the in-cloud variance of cloud water directly based on the Var , taking the results of the two PDF components as a whole. The main difference between Equation 8 and Equation 10 may come from the difference between 1/CF and ∑ 2 =1 ∕CF . As is shown in Section 3, this calculation also does not produce reasonable ν.

Cases Setting
We set two cases with different formulas of ν in warm rain simulation in CESM2.1.0 ( Package (COSP, Bodas-Salcedo et al., 2011). COSP distributes the grid-mean cloud properties from GCM into 100 subcolumns which are similar to pixels in satellite remote sensing. It also removes the subvisible clouds with cloud optical thickness (COT) below 0.3 because the passive sensor is usually unable to detect subvisible clouds (Ackerman et al., 1998). We note that each cloudy subcolumn in COSP is filled with the same amount of liquid/ice content and COT inputted from CAM6. So COSP does not allow to account for the subgrid variability of cloud water (Bodas-Salcedo et al., 2011;Klein & Jakob, 1999). The CLUBB with Subgrid Importance Latin Hypercube Sampler (CLUBB-SILHS) has been recently developed to integrate the local microphysical formula over the CLUBB PDF to upscale it to a large grid box (Larson & Griffin, 2013). However, this approach is computationally more expensive and has not been integrated with COSP yet. It would be desirable to integrate CLUBB-SILHS with COSP to enable the subgrid treatment of cloud water in COSP in the future.

Observation Data
In Sections 3.1 and 3.2, we use 10 years (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) latest collection 6 (C6) daily mean, 1°×°1 spatial resolution Level 3 cloud retrieval product of LCF, COT, and cloud droplet effective radius (CER) from the Aqua MODIS instrument (Platnick et al., 2017). MODIS level-3 provides the joint histograms of liquid-phase COT versus CER, which show the counts of successful cloud property retrievals bounded by 13 COT bin boundaries (0-150) and 10 CER bin boundaries (4-30 μm). For each COT-CER bin, liquid water path (LWP) is estimated by the following formula: where is the density of water. We can estimate the grid-mean LWP (LWP ) and the variance of LWP (Var(LWP) ) through the following equations with joint PDF of COT and CER (COT, CER) : When comparing the subgrid variance of cloud water in CAM6 and MODIS, the aforementioned methods have two limitations. First, the "one degree" spatial resolutions in CAM6 and MODIS are not exactly the same (0.9° × 1.25° vs. 1° × 1°), which may cause inconsistency to a certain extent. The data from CAM6 are regrided to the same resolution as MODIS when making comparisons. Second, MODIS only provides LWP retrievals, but the variance of subgrid cloud water in CAM6 is derived from LWC at each vertical level. As mentioned in Zhang et al. (2019), it is challenging to get the remote sensing of cloud water vertical profiles for liquid clouds from satellite sensors (Boutle et al., 2014;Lebsock et al., 2013). In this study, we assume that the variance of cloud water is not strongly related to height and use the horizontal subgrid variation of LWP to calculate ν as an observational reference for convenience, then using Equation 4 to get EF. The vertical dependence of the horizontal variance of cloud property was recently investigated in Covert et al. (2022), but it is only a case study based on in situ observations. More studies are needed in the future to achieve a more comprehensive understanding.
The Clouds and the Earth's Radiant Energy System (CERES) provides radiation measurements  at the top of the atmosphere. We use the SWCF from CERES_SYN1deg edition 4.1 hourly data from 2006 to 2015 to compare with the model results.
For the analysis in Section 3.3 about ACI, we collected the MODIS Collection 6.1 Level 2 Aqua data in the year 2015. MODIS granule was divided into 1° × 1 °scenes using the methodology developed in Zhu et al. (2018).
For the reliability of retrieved cloud properties, we rejected the pixels with CER uncertainty at 2.1 μm > 10% and solar zenith angle >65°. We only select the data of marine single layer liquid phase clouds. LCF is defined by the fraction of pixels whose reflectance at 0.61 micron is greater than the reflectance that corresponds to the minimum detectable COT. CDNC is retrieved in the convective cores of clouds (pixels with the highest 10% of COT within the scene) which are closest to adiabatic, and is used to represent the effect of aerosol (Rosenfeld et al., 2019;Zhu et al., 2018). In order to control meteorology conditions when examining ACI, we classify the data into cloud geometrical thickness (CGT) bins (Rosenfeld et al., 2019). CGT is calculated with LWP, cloud top 6 of 20 temperature (CTT) and sea surface temperature (SST) when assuming dry adiabatic lapse rate from SST to cloud base, and moist adiabatic from cloud base to cloud top. The radiative fluxes are calculated by the delta-Eddington technique in conjunction with an exponential-sum-fit parameterization (Briegleb, 1992). The incoming solar radiation is calculated by the position and the date of each grid point at local time 1:30 p.m. The height of the cloud is fixed at 800∼900 hPa when calculating radiative fluxes. The direct effect of aerosol on radiative flux is not included. To consider the subgrid variance of cloud properties in each 1° grid point, we first classify the data of COT and CER in the subcolumns into bins. Then we calculate the mean values of COT and CER in bins and feed them to the radiation code to get the mean SWCF weighted by the number in each bin. To match with the local time of MODIS Aqua data, we select the model data with local time from 12:00 at noon to 3:00 p.m. in Section 3.3, because the outputs of the model have a 3-hr interval.

Improvements in Simulated ν and EF in CLUBB
First, we focus on the improvement of ν and EF in INCLD_ν case compared to Default case and MODIS. Figure 1 shows the horizontal distribution of ν and EF in the two cases and MODIS for warm clouds (CTT > 273K). We plot the median values to avoid the influence of extremely large ν when the variance of the cloud water approaches zero. In MODIS, ν tends to be largest in the Sc decks off western coasts of continents and it diminishes significantly toward the open ocean with broken Cu. This means cloud water is more homogeneous in the Sc regions with large cloud fraction and stable meteorological conditions, and is more inhomogeneous in the Cu regions with smaller cloud fraction and more active convection events. ν in the Default case also shows the Sc-to-Cu transition pattern (Figure 1b), but it is overall underestimated over the open ocean, except western coasts of American continents and regions near 45° where the LCF tends to be large ( Figure 1c). Most of the oceanic areas have ν below 1, representing too large subgrid cloud water variance. According to Equation 4, EF increases with decreasing ν and is sensitive to ν because of the large exponent (2.47) in the denominator. Small biases in ν can cause huge differences in EF. This causes the extremely large EF over most of the open ocean in Default case, especially in Cu regions where the peak of EF can exceed 100 (Figure 1e).
The formula of in-cloud subgrid water variance in INCLD_ν case greatly improves the magnitude of ν and EF (Figures 1c and 1f). Equation 8 only uses the in-cloud cloud water in calculating, avoiding the overestimated variance caused by the clear-sky part of PDF with zero cloud water mixing ratio in a grid point. The presence of 7 of 20 LWC instead of LWCgrid in the numerator in Equation 9 also makes ν larger. This improvement is more obvious in the areas with small cloud fraction and more active convection events, like the Cu regions over the open ocean.
The new calculation of ν can also capture the observed pattern that the cloud water variance and EF are smaller in the Sc areas near the coasts. The EF over the open ocean is slightly underestimated compared to the observation.
We also test the results with Equation 10 derived in H. Guo et al. (2014). Equation 10 cannot reproduce the Sc-to-Cu transition pattern in the observation ( Figure S1 in Supporting Information S1), showing larger EF in Sc regions. The derivation of Equation 10 does not consider that the cloud properties in CLUBB come from the combination of two PDFs. Deriving the cloud properties in each PDF separately as in Equation 8 matches better with the design of CLUBB. In the rest of this study, we will focus on the improvements in microphysical processes, cloud properties, and ACI due to the updated formula in INCLD_ν.

Microphysical Processes Rates
Adjusting the EF directly changes the magnitude of the autoconversion rate. Figure 2 shows the horizontal distribution of the annual-mean vertical integrated in-cloud process rates (unit: g/m 2 /s) in liquid phase clouds. An obvious decrease in autoconversion can be found over most of the oceanic regions in the INCLD_ν case ( Figure 2b). The pattern of this change matches well with the change of EF. At the same time, the accretion rate noticeably increases, especially in Cu areas (Figure 2c). The accretion rate shows the collection rate of cloud drops by rain. It increases with larger cloud water and rainwater content. The accretion rate with prognostic precipitation in MG2 is stronger than that with diagnostic precipitation, and the accretion to autoconversion ratio increases rapidly with an increase in cloud water . The decrease in autoconversion rate leads to less consumption of cloud water. The accumulated cloud water results in the increase in accretion rate. As a result, the total conversion rate of cloud water (autoconversion + accretion) increases over Cu areas in the open ocean and decreases over Sc regions near continents and cloudy regions near northern and southern 30° (Figure 2d). Figure 3 shows the mean in-cloud autoconversion rate (unit: kg/kg/s) and EF from the model in given LWC and CDNC bins. LWC and CDNC are in-cloud values that are fed into the cloud microphysical scheme, and so they are directly related to the calculation of autoconversion rate. More cloud water and larger droplet size lead to stronger collision-coalescence of cloud droplets, meaning that autoconversion rate should be larger when LWC is larger and CDNC is smaller. However, the autoconversion rate in the Default case peaks when LWC and CDNC are both small (lower left part of Figure 3a). This unphysical phenomenon is caused by EF. According to previous studies (Hotta et al., 2020;Zhang et al., 2019), EF tends to be large in clouds with small LCF, which usually have small LWP and CDNC. As is shown in Figure 3c, EF in the Default case is overestimated most severely when LWP and CDNC are small, resulting in a distinct peak of autoconversion rate when in-cloud liquid water is low (Figure 3a). In the INCLD_ν case, the magnitude of EF decreases significantly and the contrast of EF between low LWC and high LWC is smaller (Figure 3b). As a result, the relationship between autoconversion rate, LWC, and CDNC is significantly improved. The magnitude of autoconversion rate overall decreases and peaks at the upper left part of the picture (Figure 3b).

Cloud Properties
The overly frequent rain production in GCMs is a well-known issue reported by previous studies (Franklin et al., 2013;Song et al., 2018;Stephens et al., 2010;Suzuki et al., 2015). To evaluate whether there is an improvement in the rain frequency, the probability of precipitation (POP) is shown in Figure 4. We define the POP by the fraction of liquid phase clouds that have surface rain rates larger than 0.14 mm/day from 3 hourly instantaneous data (Bai et al., 2018;Comstock et al., 2004). The POP in the Default case exceeds 70% over most of the oceanic areas, showing especially large values in Cu regions with the most severely overestimated EFs (Figure 4a). The pattern of the change in POP in INCLD_ν case is similar to that of the total conversion rate (Figure 2d). POP noticeably decreases in Sc regions near continents (Figure 4b). The relative change ((Default-INCLD_ν)/Default) is less than 5% in most of the areas where POP increases. Figure 5 shows the comparison of annual mean in-cloud COT and LCF between MODIS and CAM6 with COSP in two cases. The Default case significantly overestimates COT especially near coastal areas in low latitudes and 8 of 20 the East China Sea (Figure 5d), which is a common problem in GCMs (Song et al., 2018;Webb et al., 2001;Zhang, 2005). In the INCLD_ν case, COT increases in Sc areas (Figure 5c) caused by the decrease in autoconversion rate and the increase of LWP ( Figure S2 in Supporting Information S1). COT decreases over the Cu regions. The patten of the change in LWP is similar to that in COT and autoconversion rate. The increase of LWP in INCLD_ν make it in better agreement with the MODIS ( Figure S2 in Supporting Information S1). LCF is underestimated in typical Sc regions near the continents in the Default case ( Figure 5i) compared to MODIS. Increases of LCF in the INCLD_ν case occur in cloudy areas near 30°S and Sc regions a little away from the coast of American and African continents (Figure 5h), where cloud water increases noticeably. The changes in in-cloud COT, and LCF affect the representation of SWCF, which is highly related to the global radiative budget and climate. Figure 6 shows the comparison of annual mean SWCF between the CERES satellite retrievals and the two cases. SWCF is calculated by the difference between all-sky and clear-sky shortwave net flux at the top of the atmosphere, which is affected by cloud fraction and COT. The Default case qualitatively captures the main features and patterns of observed SWCF (Figures 6a and 6b), but still has noticeable biases ( Figure 6d). The Default case significantly underestimates the cooling effect in coastal Sc regions especially in the eastern Pacific and Atlantic, and overestimates it in other regions. In the INCLD_ν case, caused by the increase of LCF together with the changes in COT (Figures 5c and 5h), the SWCF is in better agreement with CERES in part of the oceanic areas (Figures 6c and 6d). The annual mean values of the SWCF over oceanic areas in two cases do not show a noticeable difference.
The changes of the cloud properties and SWCF over the areas poleward of 45° are not as obvious as that over lower latitudes (Figures 5c, 5h, and 6c), but SWCF is also improved in INCLD_ν in some areas. The LWP and COT increase near the north pole in INCLD_ν case, becoming closer to that in MODIS. Combined with the increase (decrease) in LCF over eastern (western) hemisphere near north pole, the SWCF matches better with CERES poleward of 60°N (Figure 6c). The SWCF over southern pacific near 60°S also becomes closer to CERES caused by the decrease of COT and LCF here. To confirm whether the change in ν and EF effects other climatologically important quantities, we calculate the annual global mean cloud and radiation parameters from 5-year model cases (without COSP) and observation whenever available (Table 2). Overall, these parameters in INCLD_ν are similar to those in the Default case.

The Improvements and Remaining Problems in Simulated ACI
The change in EF can have large effects on ACI, so here we examine the ACI in CAM6 in two cases. The study regions are the oceanic areas between 45°S and 45°N. We focus on the marine single layer liquid phase clouds that are not obscured by the higher cloud layers. Both model and MODIS data have a restriction of LCF >0.02.
The MODIS data in this part have been introduced in Section 2.4. In CAM6, we use 3 hourly instantaneous data in the years 2001∼2005. We use the maximum CDNC within vertical levels in each cloud scene to match with CDNC from the convective cores in MODIS. CGTs in models are calculated by the difference between the height of vertical interfaces at cloud base and cloud top.
In our results, over 90% of the CGTs in MODIS data are below 800 m. However, we will lose a large part of model data if we also set the maximum CGT to 800 m in two model cases. Actually, the CGTs in MODIS calculated through the assumption of adiabaticity may be underestimated because the adiabatic fraction of the cloud water is often much smaller in reality than the assumed value (Efraim et al., 2020;Lu et al., 2021). As a result, the model data with thicker clouds is comparable to MODIS when we focus on the overall comparison. To include most of the data both in model and observation, we choose the data with CGTs from 150 to 1,500 m in CAM6 and the data with CGTs from 0 to 800 m in MODIS. The data in the same CGT bin is not comparable between the model and MODIS. Stricter comparison of CGTs between MODIS and models still needs further work in the future. The mean sensitivities of cloud properties to aerosol, weighted by the amount of data in each CGT bin, are listed in Table 3.
The sensitivities of the in-cloud autoconversion rate and total conversion rate are shown in Figure 7. To match with the selection of CDNC, these process rates are the values at the vertical levels with maximum CDNC in each cloud scene in the model (unit: kg/kg/s). In Default case, the autoconversion rate in smaller CGT bins can be very large caused by the overlarge EF when LWC is small (Figure 7a). In Section 3.2, we have found that in the INCLD_ν case the EF contrast under low CDNC and high CDNC is smaller compared to that in the Default  12 of 20 case (Figures 4b and 4d). This further leads to more obvious decrease in the autoconversion rate under low CDNC and a smaller contrast under low CDNC and high CDNC. As a result, the mean sensitivity of the autoconversion rate to CDNC (S auto = dln (autoconversion rate)/dln (CDNC)) in the INCLD_ν case (Figure 7b, −0.79) is smaller than that in the Default case (Figure 7a, −1.25). The sensitivities of the total conversion rate to CDNC (S con = dln (conversion rate)/dln (CDNC)) are similar to S auto in both cases (Figures 7c and 7d), regardless of the changes in accretion rate. This is because accretion has little dependence on CDNC in contrast to autoconversion (Posselt & Lohmann, 2009), and the change in S con is controlled by the change in the autoconversion rate.
The decreases of S auto and S con weaken the sensitivity of the consumption of cloud water to aerosol. This means the suppression effect of precipitation caused by increased aerosol is weakened, which further weakens the increase in cloud water. Also, the increase in accretion and the decrease in autoconversion shift the importance of warm rain production to the former, which also weakens the ACI because autoconversion rate is more related to CDNC Posselt & Lohmann, 2009;Wang et al., 2012). As a result, there are increases in the magnitudes of cloud properties under low CDNC, and the sensitivities of cloud properties to increasing aerosols decrease in INCLD_ν case.
The sensitivities of LCF and LWP are shown in Figure 8. The mean sensitivity of LCF to aerosol (S LCF = dln (LCF)/dln (CDNC)) in CGT bins is underestimated in the Default case (0.21) compared to the value in MODIS (0.28). The mean S LCF in the INCLD_ν case (0.08) decreases because of the increase of LCF when CDNC is small (Figures 8a-8c). The obvious increase of LCF with increasing aerosols is found to have a significant contribution to cloud radiative effect in MODIS (Rosenfeld et al., 2019). However, the LCF in both cases is hard to increase when CDNC is larger than 30 cm −3 and CGT is larger than 450 m. When CGT is below 450 m, the models produce similar S LCF to MODIS, indicating that the biases mainly reside in clouds with more active convection events.
The decrease in S auto in INCLD_ν significantly weakens the response of LWP to aerosol. The mean sensitivity of in-cloud LWP (S LWP = dln (LWP)/ dln (CDNC)) is near zero in MODIS (Figure 8f, −0.04). Compared to the observation, the mean S LWP is significantly overestimated in the Default case (Figure 8d, 0.20), which is a common problem in GCMs as discussed in the introduction. The mean S LWP in INCLD_ν decreases significantly (Figure 8e, 0.07), which is in better agreement with MODIS. According to the previous research, the single-column simulations with CLUBB can produce both positive and negative LWP responses to increasing aerosol concentrations as in MODIS, depending on precipitation and free atmosphere relative humidity (Ackerman et al., 2004;H. Guo et al., 2011). With the improved EF and decreased S auto , we expect that the negative LWP responses can be shown in conditions with low relative humidity over the inversion height in CAM6.
The weakened S LWP in the INCLD_ν case affects the sensitivities of in-cloud COT and cloud top CER (Figure 9). Compared to the Default case, the CER in the INCLD_ν case increases under low CDNC (Figures 9a and 9b) and is closer to the value in observation (Figure 9c), which is associated with the increase of LWP in the fixed CDNC bins (Figure 8e). The mean sensitivity of CER (S CER = dln (CER)/dln (CDNC)) does not change much. Larger LWP is associated with larger COT, so COT also increases when CDNC is small in INCLD_ν (Figures 9d and 9e). As COT ≈ 3LWC 2 CER (Hartmann, 2015), the mean sensitivity of COT (S COT = dln (COT)/dln (CDNC)) is similar to the sum of S LWP and S CER . S COT decreases from 0.39 in Default to 0.30 in INCLD_ν caused by the decrease in S LWP . The S CER and S COT in INCLD_ν are both in better agreement with that in MODIS. The S COT in INCLD_ν and MODIS are both near 0.3, which matches well with the relationship "COT ∝ CDNC 1/3 " (Platnick & Twomey, 1994).
The effects of cloud properties on SWCF can be divided into two parts: the effect of LCF and the effect of COT. LWP and CER contribute to COT, and COT contributes to in-cloud albedo and in-cloud SWCF. In-cloud SWCF and LCF together contribute to the grid-mean SWCF. We do similar analyses on in-cloud SWCF and grid-mean SWCF ( Figure 10). The in-cloud SWCF is calculated by dividing the grid-mean SWCF by the LCF.
Both cases underestimate the sensitivity of in-cloud SWCF (S in_SWCF = dln (in-cloud SWCF)/dln (CDNC),0.25,0.19,and 0.29 in Default,INCLD_ν,and MODIS). The S in_SWCF is smaller in INCLD_ν than in Default due to the decrease in S COT . Although S COT is similar to S in_SWCF in MODIS, the S COT is smaller than S in_SWCF in two model cases. Together with the underestimation of S LCF in model, the mean sensitivities of grid-mean SWCF (S grid_SWCF = dln (grid-mean SWCF)/dln (CDNC)) are also underestimated (0.44,0.31,and 0.58 in Default,INCLD_ν,and MODIS).
The difference between S COT and S in_SWCF in CAM6 is partly caused by the overestimation in COT. For conservative scattering, which in the absence of significant aerosol absorption, the two-stream approximation used in the radiation calculation gives ln( ) ln(COT) = 2 (2+(1− )COT) and ln( ) ln(CDNC) = 2 3(2+(1− )COT) (A is the cloud albedo, g is the asymmetry parameter, Painemal, 2018;Platnick & Twomey, 1994). g is approximately constant for conservative scattering. This means the value of cloud albedo is hard to increase when COT is large, and the sensitivity of cloud albedo to CDNC decreases with increasing COT. In CAM6, the magnitude of COT is overestimated and increases fast with   . c CLDTOT is obtained from ISCCP for the years 1983-2001 (Rossow & Schiffer, 1999), MODIS data for the years 2001(Platnick et al., 2003, and HIRS data for the years 1979-2001 (Wylie et al., 2005). d LWP is derived from the Special Sensor Microwave/Imager (SSM/I, Ferraro et al., 1996;Greenwald et al., 1993;Weng & Grody, 1994) and ISCCP for the year 1987 (Han et al., 1994). SSM/I data are restricted to oceans. e PRECT is taken from the Global Precipitation Climatology Project (GPCP) for the years 1979-2003 (Adler et al., 2003, http://www.gewex.org/gpcpdata). f Radiative fluxes are from the CERES data introduced in Section 2.4. g The uncertainties for the radiative fluxes are taken from Stephens et al. (2012). No uncertainty is provided for FSNTOAC.  (Figures 9d and 9e), resulting in small sensitivity of cloud albedo and S in_SWCF to CDNC. In contrast, the COT in MODIS is generally smaller compared to that in CAM6.
This underestimation of the sensitivity of cloud albedo to CDNC is also a problem in the regional model. Liu et al. (2020) did a similar analysis on the sensitivities of cloud properties over the southeast Pacific using the Weather Research and Forecasting Model with chemistry (WRF-Chem) with the Morrison microphysical scheme (Morrison et al., 2005(Morrison et al., , 2009. The sum of the magnitude of S CER (0.16) and S LWP (0.26) is also much larger than the sensitivity of cloud albedo (0.17) in their results, indicating the necessity of further work on understanding the sensitivities of cloud properties and their radiative forcing in both regional and global models.

Summary and Discussion
This study explores the appropriate coupling between nonlinear microphysic processes and macrophysic scheme with subgrid parameterization in GCMs and its importance for simulating clouds, radiation fluxes and ACIs. We develop a new formula for the in-cloud subgrid cloud water variance in CLUBB to improve the inverse relative variance of cloud water, ν, and EF of autoconversion rate (EF) in warm rain simulation. Using MODIS as a benchmark, the new formula greatly improves the magnitude of ν and EF as compared to the case with the default formula. The new formula better captures the observed Sc-to-Cu transition. This contributes to a better simulation of cloud properties and SWCF.
With the improved EF, the S LWP decreases noticeably and is in better agreement with that in MODIS. This is caused by the decrease in S auto due to the smaller contrast in the autoconversion rate under low CDNC and 16 of 20 high CDNC. According to the single-column results (Guo et al., 2011), we expect a negative LWP response in conditions with low relative humidity in CAM6. Although the sensitivity of COT to aerosol is well produced, CAM6 underestimates the sensitivity of grid-mean SWCF due to its underestimation of the sensitivities of LCF and in-cloud SWCF. The LCF in CAM6 is hard to increase with aerosol when CGT is large. The sensitivity of in-cloud SWCF is smaller than that of COT, which is partly caused by the overestimation in the magnitude of COT. GCMs tend to overestimate the brightness of clouds to obtain reasonable cloud radiative forcing, especially when LCF and CDNC are small (Song et al., 2018;Webb et al., 2001;Zhang, 2005). This overestimation makes the cloud albedo hard to increase with the growth of cloud caused by aerosol. Our results indicate the importance of simulating reasonable LCF and COT in constraining the ACI. If the S LCF and S in_SWCF in models can be made to agree better with observations, we expect an increase in the sensitivity of cloud radiative effect to aerosol. Our results indicate the importance of the appropriate coupling between nonlinear microphysic processes and macrophysic scheme with subgrid parameterization in GCMs for climate simulations. The insufficient representation of cloud variability might partly account for the biases commonly found in GCMs in cloud properties and ACI (e.g., Lebsock et al., 2013). To better simulate the effect of subgrid cloud variance on microphysical process rates, many other factors remain to be taken into account. The different PDF assumptions in CLUBB and MG2 can contribute to the bias in EF. MG2 assumes a gamma distribution of cloud water and derives EF by gamma function. However, MG2 directly uses the ν in CLUBB, which is derived by double-Gaussian joint PDFs. Furthermore, the subgrid variance of other cloud properties can also affect microphysical process rates. The subgrid variance of CDNC in polluted regions is found to have a comparable or even larger effect on autoconversion rate than that of subgrid cloud water variance (Zhang et al., 2019). Also, global models that do not take subgrid-scale covariability of cloud and precipitation water into account in microphysics parameterizations may significantly 17 of 20 underestimate grid-mean process rates in warm clouds (Lebsock et al., 2013). As mentioned in Section 2.3, Larson and Griffin (2013) integrated the local microphysical formula over the CLUBB PDF to upscale it to a large grid box (CLUBB-SILHS), which is another way to deal with the effect of subgrid cloud variability on microphysical processes and to overcome the aforementioned problems. LWC is the grid-mean cloud water in the ith PDF; LWCgrid = 1LWC1 + 2LWC2 ; is the mean of in the ith PDF; 2 is the variance of in the ith PDF.