Reconciling Compensating Errors Between Precipitation Constraints and the Energy Budget in a Climate Model

Precipitation microphysics and the effective radiative forcing due to aerosol‐cloud interactions (ERFaci) contribute to some of the largest uncertainties in general circulation models (GCMs) and are closely interrelated. This study shows that a sophisticated, two‐moment prognostic precipitation scheme can simultaneously represent both warm rain characteristics consistent with satellite observations and a realistic ERFaci magnitude, thus reconciling compensating errors between precipitation microphysics and ERFaci that are common to many GCMs. The enhancement of accretion from prognostic precipitation and accretion‐driven buffering mechanisms in scavenging processes are found to be responsible for mitigating the compensating errors. However, single‐moment prognostic precipitation without the explicit prediction of raindrop size cannot capture observed warm rain characteristics. Results underscore the importance of using a two‐moment representation of both clouds and precipitation to realistically simulate precipitation‐driven buffering of the cloud response to aerosol perturbations.


Introduction
Aerosols play an important role in the climate system by modulating the hydrological cycle and energy budget through interactions with clouds (Albrecht, 1989). Aerosol-cloud interactions (ACIs), however, are much more complex than currently represented in general circulation models (GCMs), which rely heavily on parameterizations of cloud and precipitation processes because of their coarse spatiotemporal resolution. These limitations lead to the following two major biases in GCMs: • too light and too frequent warm rain formation (e.g., Stephens et al., 2010), and • too strong ACI (e.g., Quaas et al., 2009;Stevens, 2015).
These biases have been identified in several models (Bender et al., 2019;Ghan et al., 2016;Malavelle et al., 2017;Suzuki et al., 2015). Jing et al. (2019) explored a link between these two biases (i.e., precipitation processes and ACI magnitude) in the context of a GCM. They revealed that satellite-based tuning of precipitation processes affects the radiative forcing due to ACI in such a way that the energy balance required to reproduce historical temperature trends is violated. These coupled effects of the precipitation microphysics and energy-budget requirements suggests the presence of compensating errors in GCMs Suzuki et al., 2013). which means that precipitation is removed from the atmosphere within a single model time step (Ghan & Easter, 1992). The DIAG precipitation framework therefore generally mutes the accretion of cloud droplets by coexisting raindrops, and thus, rain formation heavily relies on the autoconversion process that directly links to aerosols. This tends to cause an exaggerated cloud water response to perturbed aerosols and thus large ACI forcing (e.g., Gettelman et al., 2015;Wang et al., 2012) as a result of the monotonic delay of precipitation onset with increasing aerosols following the autoconversion parameterization (Michibata et al., 2016). This is in stark contrast to observations that both increasing (cloud lifetime effect Albrecht, 1989) and decreasing (buffering effect Stevens & Feingold, 2009) responses of cloud water to increased aerosols can occur depending on the cloud regime Michibata et al., 2016;Toll et al., 2017Toll et al., , 2019.
To address this current limitation in GCMs, more physically based representations of aerosol-cloud-precipitation interactions are critical (Mülmenstädt & Feingold, 2018;Stephens et al., 2019). In this spirit, some GCMs have incorporated a prognostic precipitation scheme (hereafter referred to as PROG) and report improved simulations of the relative magnitudes of autoconversion and accretion (Posselt & Lohmann, 2008), precipitation rates (Abel & Boutle, 2012), hydrometeor distributions (Sant et al., 2015), and radiative forcing due to ACI .
Motivated by these pioneering works, we introduced a PROG precipitation scheme into the MIROC6 GCM (Michibata, Suzuki, Sekiguchi, et al., 2019) and obtained results that generally support the findings described above. This major update in the treatment of precipitation permits an investigation of how an improved representation of precipitation can help mitigate the aforementioned dichotomy between satellite-based constraints on warm rain process and the energy balance requirement on ACI. Here, we compare warm rain characteristics simulated using different treatments of precipitation (i.e., DIAG vs. PROG) against satellite observations by applying process-oriented metrics. We also investigate how these differences influence the cloud response to aerosol perturbations and identify links to the magnitude of the effective radiative forcing due to ACI (ERF aci ). These links are evaluated using a suite of sensitivity experiments that perturb key physical processes within a single model and are thus unaffected by intermodel differences in parameterizations.

Model Description
We use the MIROC6-SPRINTARS global aerosol-climate model (Tatebe et al., 2019) in this study. The aerosol module SPRINTARS treats the main aerosol species in the troposphere (black carbon, organic matter, sulfate, dust, and sea salt) and gas-phase precursors of sulfate (Takemura et al., 2000). The model represents cloud liquid water and ice using a two-moment bulk scheme Watanabe et al., 2009) that prognoses mass and number mixing ratios of cloud hydrometeors. Precipitation is removed from the atmosphere within a single model time step in the original version of the model that treats precipitating hydrometers diagnostically (i.e., DIAG).
We also use another version of the model, which uses a prognostic precipitation scheme (i.e., PROG) called Cloud/Hydrometeors Interactive Module with Explicit Rain and RAdiation (CHIMERRA; Michibata, Suzuki, Sekiguchi, et al., 2019). CHIMERRA prognoses rain and snow for both mass and number mixing ratios (i.e., two-moment), and also calculates the radiative effect of precipitating hydrometeors.
To examine the sensitivity of the magnitude of ACI, we test widely used autoconversion schemes summarized in Michibata and Takemura (2015) and Jing et al. (2019) using both DIAG and PROG. In PROG, an additional scheme described by Seifert and Beheng (2006) is tested, as it uses the rainwater mixing ratio provided by the PROG scheme. The results are summarized in Table S1 in the supporting information.
The model resolution is 1.4 • × 1.4 • with 40 vertical layers (T85L40). Simulations are performed for 6 years with prescribed climatological sea surface temperature and sea ice. The subsequent 5 years are used for the analysis. An evaluation of the simulated precipitation process (described in section 2.3) is carried out using the CFMIP Observational Simulator Package (COSP; Bodas-Salcedo et al., 2011;Swales et al., 2018), which enables scale-aware and definition-aware comparisons between model results and observations (Kay et al., 2018). The COSP is called every three model hours and is operated for the final year.

A-Train Satellite Data
We use the CloudSat Level 2B-TAU product (Polonsky, 2008), which is collocated with MODIS observations, the 2B-GEOPROF vertical radar product (Marchand et al., 2008), and the ECMWF-AUX reanalysis product (Partain, 2007). The analysis is restricted to single-layer warm clouds (SLWCs) for comparisons with the large-scale condensation clouds of the model. We obtained more than 7.8 million SLWCs from a full 5-year analysis from June 2006 to April 2011.
The combined use of CloudSat and MODIS data enables us to investigate cloud-to-rain conversion in the context of the vertical microphysical structure of SLWCs for different cloud regimes, which reveals modelobservation differences in warm rain characteristics (Jing et al., 2017;Suzuki et al., 2015), as described below.

Diagnostics
To illustrate how model-observation differences in rain formation process relate to the treatment of precipitation in the model, we apply metrics that probe the process in the form of the probability density function (PDF) for radar reflectivity profiles rescaled by the vertically sliced in-cloud optical depth (ICOD). This method, referred to as the contoured frequency by optical depth diagram (CFODD; Nakajima et al., 2010), is particularly useful in evaluations of how the vertical microphysical structures of warm clouds differs between nonprecipitating and precipitating regimes depending on the cloud-top effective radius R e . These diagnostics have also been implemented in COSP2 as an inline warm-rain diagnostic tool (Michibata, Suzuki, Ogura, et al., 2019), which is employed in this study.
This study applies clean-sky ERF aci to exclude contamination due to aerosol scattering and absorption in cloudy-sky conditions; it is defined as the difference in cloud radiative forcing between present-day (PD) and preindustrial (PI) aerosol emissions under the clean-sky (i.e., aerosol-free) condition (Ghan, 2013). All simulations ensure that the imbalance in the radiative flux at the top-of-the-atmosphere (TOA) is within 1.0 W m −2 under the PD condition. This is achieved by parameter tuning for warm rain efficiency alone, and a suite of sensitivity experiments was conducted in which only the warm rain parameters were varied. Figure 1 compares warm rain characteristics in the form of CFODD statistics between A-Train satellite observations and MIROC6 with the DIAG and PROG schemes. As previous studies have reported , GCMs with traditional DIAG schemes typically simulate precipitation events too frequently, even in the smallest R e range (Figure 1d) where satellite observations indicate no precipitation (Figure 1a). This is a well-known problem in most GCMs (Jing et al., 2017;Suzuki et al., 2015). However, the PROG model well captures the nonprecipitating regions (<-15 dBZe) in the smallest R e range (Figure 1g). The mode of the PDF also shifts toward higher dBZe values with increasing R e through drizzling (-15 dBZe <) to precipitating (0 dBZe <) characteristics (Figures 1h and 1i), in better agreement with satellite-derived CFODD statistics (Figures 1a-1c).

Reconciling Compensating Errors Between Microphysics and ACI
This systematic difference of CFODDs between DIAG and PROG can be interpreted as a result of the improved representation of the relative magnitudes of autoconversion and accretion. The DIAG model cannot explicitly represent the accretion of cloud droplets by coexisting raindrops because the model does not retain precipitation information across model time steps. Its rain formation thus solely relies on autoconversion for all cloud regimes (Michibata, Suzuki, Sekiguchi, et al., 2019). However, the PROG scheme can accelerate rain formation by enhancing accretion in precipitating regimes, thus reducing the relative contribution of autoconversion compared with the DIAG scheme (e.g., Gettelman et al., 2015;Posselt & Lohmann, 2008). The PROG scheme indeed shows a lower absolute value of the autoconversion rate (P aut ) than the DIAG scheme (see also Table S1) and thus mitigates the premature onset of rain, which enables the PROG model to better capture the nonprecipitating regime. This improvement in simulated rain formation is also found in geographical distributions of the fractional occurrences of warm rain regimes ( Figure S1), suggesting that the overly frequent warm rain bias (e.g., Stephens et al., 2010) is significantly reduced in the PROG model.
The CFODD statistics for the PROG model show wider PDF distributions that are closer to A-Train observations than are those of the DIAG model. Given that radar reflectivity is proportional to the sixth power of raindrop size (Stephens & Haynes, 2007), larger raindrops that form via enhanced accretion in the PROG model can result in higher dBZe values. In this context, we explore how the CFODD statistics are sensitive to the modeled raindrop size ( Figure S2). This implies that the variation in raindrop size from less than 50 m to more than 500 m, depending on the cloud regime characterized by the R e range, is required to replicate the observation-like warm rain formation process. As the PROG model includes a two-moment representation of precipitation, the size of precipitating particles is explicitly predicted, which leads to larger variability of the simulated radar reflectivity and thus a wider PDF. This is confirmed by the sensitivity of the CFODDs to whether the model prognoses raindrop number concentration (N r ) or not with the aid of a single-moment version of the PROG scheme. The single-moment PROG results in CFODD statistics that are worse than the original two-moment PROG ( Figure S3), demonstrating that the two-moment representation of precipitation is critical to simulating observed cloud-to-rain conversion.
Nevertheless, it is interesting that CFODDs with the single-moment scheme are better than the two-moment model in the large R e (>18 m) category. The single-moment scheme predicts N r significantly lower (at least an order of magnitude) than in the two-moment scheme. This is probably because the single-moment scheme does not retain N r across multiple time steps, resulting in larger raindrops characterized by higher dBZe. This does not contradict the previous studies based on the single-moment prognostic precipitation framework (e.g., Abel & Boutle, 2012). This suggests that the two-moment scheme underestimates the observed upper bound of raindrop sizes. Future studies should consider the development of more accurate raindrop size distributions (e.g., Paukert et al., 2019).
To explore how rain formation process affects macroscopic energy budget perturbations due to aerosols, we quantify the magnitude of radiative forcing due to ACI for DIAG and PROG (Figure 2). In Figure  Figure 2. Global annual mean ERF aci estimated from observation-based studies (black; Chen et al., 2014;Christensen et al., 2016Christensen et al., , 2017Douglas & L'Ecuyer, 2020) and their probable range calculated by correcting the effect of retrieval limitations (box-whisker) and that simulated by MIROC6 with various precipitation and autoconversion schemes (colors). Symbols and error bars for MIROC6 represent the mean and standard deviation of the interannual variability, respectively. The shaded area represents the uncertainty range of ERF aci estimated from IPCC AR5 (Boucher et al., 2013;Myhre et al., 2013), and the horizontal dashed line indicates the best estimate of IPCC AR5 (-0.45 W m −2 ). 2, although ERF aci values from multisatellite studies (Chen et al., 2014;Christensen et al., 2016Christensen et al., , 2017Douglas & L'Ecuyer, 2020) are also shown for reference, the observations probably underestimate the magnitude of ACI due to the inherent uncertainty in retrieval limitations (Ma et al., 2018). We therefore also provide a new estimate of multisatellite ERF aci incorporating retrieval corrections (see Text S1 in the supporting information),which shows more negative ERF aci . This is consistent with a recent satellite-based study (Rosenfeld et al., 2019) that shows stronger ACI than previous estimates from satellites. Nevertheless, the DIAG results have a more negative ERF aci than satellite-based estimates that varies with the autoconversion scheme used in the model. Results from two autoconversion schemes-those of Khairoutdinov and Kogan (2000) and Beheng (1994), which are both characterized by a strong nonlinear dependence on cloud droplet number concentration (N c )-lie below the lower bound of the uncertainty range of the IPCC AR5 (Boucher et al., 2013;Myhre et al., 2013), although these schemes tend to perform better in terms of warm rain formation . This is because rain formation process is fundamentally related to ACI via their mutual dependence on N c Wood, 2005). The global annual mean ERF aci from four autoconversion schemes in the DIAG model is -1.42 W m −2 , which is more than three times the best estimate of the IPCC AR5 (-0.45 W m −2 ; Boucher et al., 2013;Myhre et al., 2013) and outside the lower bound of multisatellite estimates (-1.38 W m −2 ).
On the other hand, the ERF aci values simulated by the PROG model are lower than those of DIAG and fall within the range estimated from satellite measurements. The global annual mean ERF aci from five autoconversion schemes in the PROG model is -0.80 W m −2 , which is within the range that can represent the observed historical temperature trend (Bellouin et al., 2020;Stevens, 2015). This means that the PROG model can simultaneously satisfy both the satellite-based constraint on precipitation microphysics and the energy budget requirements. The PROG results appear to be insensitive to the autoconversion scheme used, suggesting the presence of some mechanism that stabilizes the simulated ERF aci through weakening cloud susceptibility to aerosol perturbations, such as the buffering effects (Stevens & Feingold, 2009) discussed below.

ACI Buffering Mechanisms
A reduction in ACI due to the use of a PROG scheme has been reported for some GCMs (e.g., Gettelman et al., 2015;Posselt & Lohmann, 2008). In MIROC6, the PROG scheme simulates the ratio of accretion to autoconversion (13.6) as approximately an order of magnitude or more higher than that for the DIAG model (0.14), consistent with previous studies (Gettelman, 2015;Michibata, Suzuki, Sekiguchi, et al., 2019). The weaker ACI in the PROG model can therefore be attributed to an enhancement of accretion relative to autoconversion because autoconversion is the only pathway through which aerosols can affect rain formation (Wood, 2005). We next investigate in-depth process-level reasons for the weakening of ACI in PROG to understand the mechanisms underlying improvements in the compensating errors described above.
Given that PROG retains rainwater information across model time steps, in contrast to DIAG, the accretion process in PROG directly responds to changes in cloud water content because the accretion rate (P acc ) depends on the product of cloud water (q c ) and rainwater (q r ) mass mixing ratios. Figure 3 illustrates how the P acc responds differently to changes in q c using various precipitation schemes. The change in accretion in PROG is more correlated with changes in q c (Figure 3b) than that in DIAG (Figure 3a). This means that an increase in q c in response to increasing aerosols due to the cloud lifetime effect is partly buffered by the enhancement of accretion in PROG, which is not present in DIAG (Figure 4).
To quantify the impact of this mechanism on ACI, we conduct an additional sensitivity experiment using the PROG model that prevents P acc from exceeding P aut ; that is, P acc /P aut < 1, referred to as the accretion-limited (ACCLMT) experiment (Figure 3c). The ACCLMT experiment shows a weak correlation between the change in q c and that in P acc , similar to the DIAG model (Figure 3a), with a magnitude of ACI that is stronger (-1.10 W m −2 ; blue open diamond in Figure 2) than that of PROG control version (PROG CTRL; -0.79 W m −2 ). The CFODD statistics are also worse than those of the PROG CTRL, in particular for the larger R e Figure 4. Schematic of the modeled cloud system response to aerosol perturbations using (a) the DIAG scheme and (b) the PROG scheme. Thick red and blue allows between processes indicate positive and negative feedbacks to the initial aerosol perturbations, respectively. The subscripts for mass (q) and number (N) mixing ratios, that is, "a," "c," and "r," mean aerosol, cloud, and rain, respectively. The P aut , P acc , and P wscav are process tendencies with regard to autoconversion, accretion, and wet-scavenging, respectively. category ( Figure S4c), in which dBZe is underestimated compared with the satellite statistics. This is because raindrops cannot grow larger due to the inhibition of accretion, suggesting that a more physically based representation of precipitation leads to better estimates of ERF aci .
Accretion also decreases N c via depletion of cloud droplets by raindrops, the so-called coalescence scavenging (Wood, 2006), which can also buffer (i.e., a negative feedback) the cloud susceptibility to aerosols McCoy et al., 2020) as schematically illustrated in Figure 4. To evaluate the impact of this process on ACI, we also conduct another experiment that fixes the tendency of N c due to accretion at zero (dN c | acc =0). This simulation also results in a more negative ERF aci (-1.04 W m −2 ; blue closed diamond in Figure 2) than does PROG CTRL (-0.79 W m −2 ), suggesting that the magnitude of the buffering of ERF aci due to this accretion-driven negative feedback reaches 0.25 W m −2 . This explains nearly 40% of the ERF aci difference between DIAG (-1.42 W m −2 ) and PROG (-0.80 W m −2 ), underscoring the significant role of coalescence scavenging in climate simulations, which is captured only when using prognostic precipitation. Although this paper focuses only on warm rain processes, ice and snow microphysical processes may also contribute to buffering in ERF aci (Michibata et al., 2020), which explain the remaining part of the ERF aci difference.
The enhancement of accretion in response to aerosols also affects the wet-scavenging efficiency, which is a key process in aerosolcloudprecipitation interactions (Figure 4). Autoconversion parameterization decreases P aut with increasing aerosol loading, which results in less efficient wet scavenging and leaves more aerosols suspended in the atmosphere, leading to a stronger ERF aci (Jing & Suzuki, 2018). This positive feedback on ERF aci from wet scavenging represented in DIAG (Figure 4a) is significantly reduced by the enhancement of accretion in response to the increased q c in PROG (Figure 3b), which leads to a negative feedback on ERF aci from wet-scavenging (Figure 4b). This is confirmed by another sensitivity experiment that fixes q c only for accretion used in the wet-scavenging parameterization. The simulation without buffering effects from wet-scavenging indeed results in a larger ERF aci (-0.92 W m −2 ) than does PROG CTRL (-0.79 W m −2 ), but the magnitude of this buffering is smaller than that due to in-cloud coalescence scavenging (-1.04 W m −2 ; see Figure 2). This suggests that coalescence scavenging of cloud droplets is more significant than wet-scavenging of aerosol particles in the buffering of ACI in MIROC6, possibly because the former directly reduces N c , which affects autoconversion, whereas the latter primarily affects the aerosol number concentration.
The results described above demonstrate that the accretion-driven buffering of ACI in PROG makes a large contribution to reducing the compensating errors related to precipitation and the energy budget.

Summary and Conclusions
This study explored the fundamental mechanism responsible for compensating errors between precipitation microphysics and energy budget requirements on ACI (Jing & Suzuki, 2018) by comparing two different treatments of precipitation, DIAG and PROG. The latter is found to mitigate the compensating errors. Specifically, we found that the PROG model can simulate cloud-to-rain transition in good agreement with satellite statistics (Figure 1) except for the larger R e (>18 m) category. In PROG, this process-based constraint on warm rain process preserves an ERF aci magnitude that is consistent with satellite-inferred estimates (Chen et al., 2014;Christensen et al., 2016Christensen et al., , 2017Douglas & L'Ecuyer, 2020) without increasing the magnitude of ACI (Figure 2). This is in a stark contrast to the more traditional DIAG model, which cannot realistically simulate precipitation and ERF aci magnitude simultaneously .
The realistic rain formation in PROG can be attributed to the improved autoconversion-to-accretion ratio for various cloud regimes Posselt & Lohmann, 2008). Incorporating a two-moment representation of precipitating particles allows the model to represent explicit changes in raindrop size, leading to variability of the simulated radar reflectivity ( Figure S2). Indeed, an investigation of the sensitivity of warm rain to whether N r is prognostic (i.e., two-moment) or not (i.e., single-moment) showed that the CFODD statistics of a single-moment version of PROG are worse than those of the two-moment version ( Figure S3). Nevertheless, the two-moment PROG underestimated the frequency of higher dBZe in a more precipitating regime (R e >18 m), implying smaller raindrop sizes than in the observations. This suggests that a more physically based representation of raindrop size distributions (e.g., Paukert et al., 2019) is critical to modeling realistic rain formation. The significant weakening of ERF aci in PROG arises from an enhancement of accretion due to the coexistence of q c and q r , which is not captured by DIAG (Figure 3). An accretion-driven buffering effect, particularly related to coalescence scavenging (Wood, 2006(Wood, , 2012, acts as a negative feedback (Figure 4) that buffers the cloud water response to aerosol perturbations. This mechanism is a major driver of the differences in ERF aci between the two precipitation schemes. The simulated ERF aci from the single-moment (-0.98 W m −2 ) and two-moment (-0.79 W m −2 ) schemes were almost the same because the enhancement of accretion that buffers the ACI is basically common to both schemes.
Results presented here underscore the importance of using a two-moment representation of both clouds and precipitation to realistically simulate rain formation and precipitation-driven feedbacks in clouds, in agreement with previous process modeling studies (e.g., Yamaguchi et al., 2017). Given that only a few GCMs incorporate two-moment prognostic precipitation schemes (Gettelman et al., , 2019Michibata, Suzuki, Sekiguchi, et al., 2019;Sant et al., 2015), future studies should investigate the sensitivity of the impact of precipitation on the magnitude of ACI found in this study to model architecture as more GCMs incorporate prognostic precipitation schemes. This paper focused on contributions from rapid adjustment to ACI, which is more likely influenced by precipitation modeling than is the Twomey effect. However, it is important in future studies to decompose the ACI into the two effects  for better quantifying the impact of precipitation modeling on relative contribution from the two effects, which is an essential source of intermodel spread of ACI .

Data Availability Statement
The CloudSat data products can be obtained from the CloudSat Data Processing Center at CIRA/Colorado State University (https://www.cloudsat.cira.colostate.edu). The A-Train data and raw simulation outputs of the MIROC6-SPRINTARS used in this study are available at https://doi.org/10.5281/zenodo.3596222 website.