Features of mid- and high-latitude low-level clouds and their relation to strong aerosol effects in the Energy Exascale Earth System Model version 2 (E3SMv2)

The E3SMv2 model, like various other global models that include representations of aerosol–cloud interactions, uses an empirically chosen lower bound on the simulated in-cloud cloud droplet number concentration (CDNC) to help constrain the effective radiative forcing of anthropogenic aerosols, ERF_aer. This study identifies where ultra-low CDNCs (i.e., concentrations lower than 10 cm⁻³) occur in the stratiform and shallow convective clouds simulated by E3SMv2 and which of the occurrences have the strongest impact on ERF_aer. Process-level analyses are presented to reveal characteristics of the cloud droplet formation and removal processes associated with impactful ultra-low CDNCs.

Simulations performed with present-day emissions show that ultra-low CDNCs are most frequently found over the mid- and high-latitude oceans in both hemispheres, while the occurrences are also frequent in polluted continental regions despite the high aerosol concentrations. Ultra-low CDNCs with the largest impacts on the simulated regional and global mean ERF_aer are found in the lower troposphere in the Northern Hemisphere middle and high latitudes. These cases are typically associated with large cloud fractions, strong water vapor condensation, weak turbulence, and lack of cloud droplet nucleation from aerosol activation. Under such atmospheric conditions, boosting aerosol activation and enhancing turbulent mixing of cloud droplet number can increase the simulated CDNCs, although the magnitude of the global mean ERF_aer increases undesirably. The reason for this model behavior is discussed. Overall, our study suggests that mid- and high-latitude low-level stratus occurring under weak turbulence is a cloud regime worth further investigating for the purpose of identifying and addressing the root causes of ultra-low CDNCs and strong ERF_aer in E3SM.

1 Introduction

The effects of anthropogenic aerosols on the energy budget of the Earth, especially the effects involving changes in cloud properties and cloud life cycle, continue to be a major challenge in our understanding of the Earth system and our ability to numerically model it (see, e.g., summaries in Seinfeld et al., 2016; Bellouin et al., 2020). The Energy Exascale Earth System Model (E3SM; Leung et al., 2020, https://e3sm.org/, last access: 31 August 2025) is a global numerical model developed for research and predictions related to Earth system variabilities. In the publicly released E3SMv1 and v2 (Golaz et al., 2019, 2022), the global mean effective radiative forcing of anthropogenic aerosols, denoted hereafter by ERF_aer, was −1.65 and −1.52 W m⁻², respectively, when estimated using pre-industrial simulations (Golaz et al., 2022, see the end of Sect. 4.2 therein). Both values were deemed large in magnitude, compared with the results from various other models (e.g., Smith et al., 2020; Zhang et al., 2022a), as well as estimates based on satellite retrievals (e.g., Bellouin et al., 2020).

Zhang et al. (2022a) showed that the strong ERF_aer in E3SMv1 was attributable primarily to the aerosol indirect effect. Golaz et al. (2022) noted that, similar to the results reported for various other models, e.g., Lohmann et al. (2000), Hoose et al. (2009), Wang and Penner (2009), and Neubauer et al. (2019), the magnitude of the aerosol indirect effect (and hence ERF_aer) in E3SM could be pragmatically reduced by imposing a lower bound on the in-cloud cloud droplet number concentration (CDNC), i.e., by letting

$$\begin{array}{}\text{(1)}& {\mathrm{CDNC}}_{\mathrm{b}}=max(\mathrm{CDNC},{\mathrm{CDNC}}_{\min }),\end{array}$$

where CDNC is the in-cloud droplet number concentration obtained by solving the model's prognostic and diagnostic equations, CDNC_min is an empirically chosen constant of the order of 10 to a few tens of droplets per cubic centimeter of air, and CDNC_b is the bounded value used in the calculation of cloud microphysics. (Since the atmosphere component of E3SM allows for partial cloudiness in a grid box – as is commonly done in global atmospheric models – we clarify that all CDNCs mentioned in this paper are the in-cloud values in the respective grid boxes. Whether the other physical quantities discussed in the following are in-cloud or grid box averages is clarified when the quantities are introduced.)

On the one hand, low CDNCs have been observed in the real world and examined by researchers: see the examples listed in Hoose et al. (2009) and more recent evidence from, e.g., Wood et al. (2018), O et al. (2018), and Grosvenor et al. (2018). On the other hand, lower bounds of CDNC were not only used in earlier simulations of aerosol indirect effects (see examples summarized in Hoose et al., 2009) but continue to be used in most recent models, e.g., 40 cm⁻³ in Salzmann et al. (2022), 20 cm⁻³ in van Noije et al. (2021), 10–20 cm⁻³ in Mignot et al. (2021), and 10 or 40 cm⁻³ in Neubauer et al. (2019). From a model developer's perspective, in order to understand the strong aerosol indirect effect in E3SM from its roots, it will be useful to understand where and why very low CDNCs occur in the simulations. To that end, the present paper combines sensitivity experiments with online and offline diagnostics to answer the “where” question about the occurrences of ultra-low CDNCs (i.e., concentrations lower than 10 cm⁻³) in E3SMv2 and provide some initial clues to the “why” question. Specifically, we identify the geographical locations and altitudes with frequent occurrences of ultra-low CDNCs and determine which of those locations have the largest impacts on the simulated ERF_aer. The features of the atmospheric conditions, as well as cloud droplet formation, transport, and removal processes associated with ultra-low CDNCs, are also discussed.

The remainder of the paper is organized as follows. Section 2 provides a description of E3SMv2's atmosphere component, focusing on the cloud and aerosol process parameterizations most relevant to this study. Section 3 presents an overview of the numerical experiments. Section 4 identifies where ultra-low CDNCs occur in E3SMv2 simulations and which of the occurrences have the strongest impact on ERF_aer. Section 5 presents process-level analyses of the atmospheric conditions and cloud droplet formation, transport, and removal processes associated with the simulated ultra-low CDNCs. Section 6 presents some additional sensitivity experiments and discussions. Section 7 summarizes the results and draws conclusions.

Figure 1 (a) Atmospheric processes considered in EAMv2 and the sequence of calculations during the model's time integration. Gray boxes correspond to processes that can directly change the mass or number mixing ratios of cloud droplets in stratiform and shallow convective clouds. (b, c) Time mean tendencies of the vertically integrated grid box mean droplet mass and number, respectively, averaged over the middle and high latitudes. The filled boxes in (b, c) show annual averages, and the whiskers indicate the ranges of monthly averages. “Cu” refers to cumulus convection. The budgets were diagnosed using 1 year of model output from a present-day EAMv2 simulation conducted without the prescribed lower bound of CDNC, i.e., configuration nc00 in Table 1.

2 Model description

The atmosphere component of E3SMv2, commonly referred to as EAMv2 (see Sect. 2.1 in Golaz et al., 2022), is a global hydrostatic atmospheric model including comprehensive representations of large- and small-scale fluid dynamics, clouds, radiation, and aerosol life cycles, as well as the interactions among these physical processes (Fig. 1a). Broadly speaking, the parameterizations used in EAMv2 divide the simulated clouds into two categories: (I) deep convection and (II) shallow convection and stratiform clouds. The two categories of cloud can co-exist in a model grid box but the latter is where the CDNCs are bounded so as to reduce the magnitude of the global mean ERF_aer and hence is the focus of this study.

To provide necessary background information for the reasoning in this paper and facilitate the comprehension of our analyses, the remainder of this section includes an overview of the sources and sinks of cloud droplets in stratiform and shallow convective clouds in EAMv2 (Sect. 2.1), a brief summary of the representation of aerosol life cycles (Sect. 2.2), and a description of the parameterization that represents the nucleation, evaporation, and turbulent mixing of cloud droplet number (Sect. 2.3).

2.1 Sources and sinks of cloud droplets

For stratiform and shallow convective clouds, EAMv2 uses the parameterization of Gettelman and Morrison (2015), commonly referred to as MG2, to represent the cloud microphysics processes, i.e., collision-coalescence of cloud condensate and formation of precipitation, as well as gravitational settling of hydrometeors. The parameterization uses a two-moment approach, in which the mass and number mixing ratios of cloud droplets are predicted with separate prognostic equations.

Among the atmospheric processes considered in EAMv2, those depicted with gray boxes in Fig. 1a can directly affect cloud droplet mass or number in stratiform and shallow convective clouds. The respective tendencies (i.e., rates of change with respect to time) are shown in Fig. 1b for droplet mass and in Fig. 1c for droplet number, presented here as the monthly and annual mean vertical integrals averaged over the middle and high latitudes. (The low latitudes, 30° S to 30° N, are excluded because the deep convective clouds dominate in those regions.) In terms of such global-scale integrals, the MG2 cloud microphysics parameterization leads to net sinks of both the mass and the number of cloud droplets (Fig. 1b and c). The primary source of cloud droplet mass is the condensation of water vapor calculated by the parameterization named CLUBB (Cloud Layers Unified by Binormals; Larson, 2017; Larson and Golaz, 2005; Golaz et al., 2002; Larson et al., 2002), which uses higher-order closures for parameterizing turbulence. The primary source of cloud droplet number is droplet nucleation (Fig. 1c), which is represented with a comprehensive parameterization following Ghan and Easter (2006) and Easter et al. (2004), as elaborated in Sect. 2.3. Condensed water detrained from deep convection is also a source for droplet mass and number for stratiform and shallow convective clouds but this source is small in the mean budget in middle and high latitudes (Fig. 1b and c). Transport by the resolved winds (i.e., large-scale advection) and by the parameterized updrafts and downdrafts in deep convection can lead to local changes of cloud droplet mass and number but the net effects after temporal and spatial averaging are very small (Fig. 1b and c).

Geographical distributions of the vertically integrated annual mean budget terms are shown in Fig. S1 in the Supplement. There, we see again that the impacts of resolved transport and parameterized deep convective transport are small. In the middle and high latitudes, the primary sources of cloud droplet mass and number are water vapor condensation and droplet nucleation, respectively, while cloud microphysics is the primary sink for both droplet mass and droplet number.

2.2 Aerosols

The abundance of aerosol particles in the atmosphere is simulated in EAMv2 by solving time evolution equations of particle mass and number mixing ratios. The wide ranges of particle size and composition are taken into account by considering seven chemical components (sulfate, black carbon, primary organic matter, secondary organic matter, marine organic matter, dust, sea salt) and using four lognormally distributed modes to represent the statistical distributions of particle size and composition at specific time instances and locations. Furthermore, aerosol particles in EAMv2 are divided into two subpopulations of different attachment states: the interstitial aerosols are those found outside cloud droplets, while the cloud-borne aerosols are those embedded in cloud droplets. EAMv2 solves separate mass and number mixing ratio equations for the different compositions, modes, and attachment states listed here. Detailed descriptions of the aerosol life cycles in EAMv2 can be found in Wang et al. (2020), Liu et al. (2016), and Liu et al. (2012).

The activation of aerosol particles leads to the formation of new cloud droplets, a process commonly referred to as cloud droplet nucleation. From the aerosol life cycle perspective, activation changes the attachment state of an aerosol particle from interstitial to cloud-borne, while the evaporation of cloud droplets can result in aerosol resuspension, i.e., the conversion of cloud-borne aerosols back to the interstitial state. Further details of the parameterization in EAMv2 involving aerosol activation can be found in Sect. 2.3.

2.3 Nucleation, evaporation, and turbulent mixing of cloud droplets

In EAMv2, the changes in cloud droplet number mixing ratio caused by cloud droplet nucleation or evaporation are represented by a scheme based on Ghan and Easter (2006) and Easter et al. (2004), similar to the summary in Sect. S1.1.8 of Liu et al. (2012). The following features of the implementation in EAMv2 (and also in EAMv1; Rasch et al., 2019) are important for understanding the results presented in later sections.

First, the nucleation, evaporation, and turbulent mixing of cloud droplets are tightly coupled both physically and numerically. For brevity, we refer to the parameterization in EAMv2 that describes this collection of processes as DROPMIXNUC, following the Fortran subroutine name in the EAMv2 code. EAMv2 does not assume vertically constant CDNCs in multi-layer cloudy regions; rather, the parameterized turbulence is responsible for the upward transport of cloud droplets formed at the cloud base and for the subgrid-scale mixing of droplets across vertical layers. This design choice is different from the assumptions used in the ECHAM4 model described in Sect. 2.3 of Lohmann et al. (1999) and in various more recent models, such as ECHAM-HAM (Zhang et al., 2012; Tegen et al., 2019; Neubauer et al., 2019), ECHAM-SALSA (Kokkola et al., 2018), and UKESM1 (Grosvenor and Carslaw, 2020). The implementation in these other models can be interpreted as assuming that turbulent transport is sufficiently strong to vertically homogenize CDNC in a contiguous cloudy region within the same time step of the numerical simulation.

Second, droplet nucleation in EAMv2 (and also v1) is treated differently at the cloud base and inside newly formed clouds, see Sect. 2.1 in Ghan and Easter (2006) as well as Sect. 2.3.2 versus Sect. 2.3.3 in the following. In both cases, however, the fractions of interstitial aerosol particles in different lognormal modes that are activated to form cloud droplets are calculated using the parameterization of Abdul-Razzak and Ghan (2000, hereafter ARG2000), assuming that the supersaturation of water vapor that leads to aerosol activation is caused by the adiabatic ascents of air parcels. In EAMv2 (and also v1), the characteristic speed of such ascents, $w_{\text{act}}^{\ast }$, is estimated as

$$\begin{array}{}\text{(2)}& {\displaystyle }{w}_{\text{act}}^{\ast }& {\displaystyle }=max(w_{\text{act,min}},{\mathit{\sigma }}_{\mathrm{w}}),\text{ where}\text{(3)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{w}}& {\displaystyle }=min\left({\mathit{\sigma }}_{\text{max}},\sqrt{\stackrel{\mathrm{‾}}{w^{\prime \mathrm{2}}}}\right).\end{array}$$

Here $\stackrel{\mathrm{‾}}{w^{\prime \mathrm{2}}}$ is the subgrid variance of the vertical velocity of air calculated by CLUBB. The minimum ascent velocity, w_act,min, is set to 0.1 m s⁻¹, and σ_max = 10 m s⁻¹.

Overall, the impact of the DROPMIXNUC parameterization on cloud droplet number can be understood as a grid box mean droplet number tendency, expressed as

$$\begin{array}{}\text{(4)}& {\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{DROPMIXNUC}}={\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{mix}}+{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{nuc-evap}},\end{array}$$

where N_l denotes the droplet number mixing ratio and an overline denotes the spatial average over a grid box. The first right-hand-side term in Eq. (4) is turbulent mixing parameterized using the eddy diffusivity theory. The second term is formulated as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{nuc-evap}}& ={\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-regen}}\\ & +{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-shrink}}+{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-grow}}\\ & +{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-base-nuc}}+{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{vert-detr}}.\end{array}\end{array}$$

Figure 2 Treatment of different scenarios and physical processes in the parameterization of cloud droplet nucleation, evaporation, and turbulent mixing in EAMv2. More detailed descriptions can be found in Sect. 2.3.

The first term on the right-hand side of Eq. (5) represents droplet nucleation by cloud regeneration inside persisting clouds (see Sect. S1.1.8 in Liu et al., 2012) but is neglected in EAMv2 (and v1). The other four terms represent droplet evaporation associated with cloud fraction decrease (Sect. 2.3.1 and Fig. 2a), droplet nucleation associated with cloud fraction increase (Sect. 2.3.2 and Fig. 2b), nucleation at the cloud base (Sect. 2.3.3 and Fig. 2c), and the evaporation of droplets detrained by turbulence (Sect. 2.3.4 and Fig. 2d). The cloud fraction used in DROPMIXNUC is the liquid cloud fraction calculated by CLUBB.

2.3.1 Droplet evaporation associated with cloud fraction decrease

If the cloud fraction in a grid box has decreased from the previous time instance at which the DROPMIXNUC parameterization was calculated, the decrement in cloud fraction is assumed to correspond to disappeared clouds (Fig. 2a). The decrement in the grid box mean droplet number mixing ratio is assumed to be proportional to (I) the change in cloud fraction and (II) the in-cloud mean droplet number concentration of the previous time instance, namely,

$$\begin{array}{}\text{(6)}& {\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-shrink}}\mathrm{\Delta }t=\widehat{N_{\mathrm{l}}}\left(t\right)\left[f_{\mathrm{c}}(t+\mathrm{\Delta }t)-f_{\mathrm{c}}\left(t\right)\right],\end{array}$$

where Δt is the time step size used by the parameterization and $\widehat{N_{\mathrm{l}}}\left(t\right)$ is the in-cloud mean droplet number mixing ratio at time instance t, while f_c(t) and f_c(t+Δt) are cloud fractions at the old and new time instances, respectively. As a result of this evaporation, the previously cloud-borne aerosol particles are resuspended, resulting in increases in the grid box mean interstitial aerosol number and mass concentration.

2.3.2 Droplet nucleation associated with cloud fraction increase

If the cloud fraction in a grid box has increased from the previous time instance, the increment is assumed to correspond to newly formed clouds (Fig. 2b). For these new clouds, the ARG2000 parameterization is used to calculate the fractions of activated aerosol particles, denoted here as f_a,i, $i=\mathrm{1},\mathrm{\dots },M$, with M being the total number of aerosol modes. Note that f_a,i depends critically on $w_{\text{act}}^{\ast }$ as well as on the aerosol properties. These fractions are applied to the mean interstitial aerosol number mixing ratios in the grid box, denoted here by ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}$, $i=\mathrm{1},\mathrm{\dots },M$. Each activated aerosol particle is assumed to correspond to one newly nucleated cloud droplet. This leads to an increase in the grid box mean droplet number mixing ratio that is proportional to the cloud fraction increment, namely,

$$\begin{array}{}\text{(7)}& {\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-grow}}\mathrm{\Delta }t=\left[f_{\mathrm{c}}(t+\mathrm{\Delta }t)-f_{\mathrm{c}}\left(t\right)\right]\sum _{i=\mathrm{1}}^M{f}_{\mathrm{a},i}\left(t\right){\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}\left(t\right).\end{array}$$

2.3.3 Droplet nucleation at cloud base

A grid box is assumed to have a fraction of its horizontal area covered by cloud base if the cloud fraction in the current grid box is larger than the cloud fraction in the grid box below by at least 0.01 (Fig. 2c). The cloud fraction difference between the two layers (a positive number) is defined as the cloud base fraction, f_b, of the grid box. At the cloud base, the fractions of aerosol particles being activated from the below-base layer are calculated using the parameterization of ARG2000. These activation fractions are applied to the turbulent fluxes of interstitial aerosols across a unit area at the cloud base, ${\stackrel{\mathrm{‾}}{\mathcal{F}}}_{\mathrm{a},i}$, where

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{‾}}{\mathcal{F}}}_{\mathrm{a},i}=w_{\text{act}}^{\ast }{\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i},i=\mathrm{1},\mathrm{\dots },M.\end{array}$$

Each activated aerosol particle is assumed to correspond to one newly nucleated cloud droplet. The resulting grid box mean droplet number tendency is

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{N_{\mathrm{l}}}}{\partial t}}\right)}_{\text{cld-base-nuc}}& ={\displaystyle \frac{f_{\text{b}}}{\mathrm{\Delta }z}}\sum _{i=\mathrm{1}}^M{f}_{\mathrm{a},i}{\stackrel{\mathrm{‾}}{\mathcal{F}}}_{\mathrm{a},i}\\ & ={\displaystyle \frac{f_{\text{b}}}{\mathrm{\Delta }z}}\sum _{i=\mathrm{1}}^M{f}_{\mathrm{a},i}{w}_{\text{act}}^{\ast }{\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i},\end{array}\end{array}$$

where z is the geopotential height and Δz is the thickness of the model layer. The aerosol number concentrations ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}$ in Eqs. (8) and (9) are the concentrations in the model layer right below the cloud base.

The droplet number increment expressed in Eq. (9) and the turbulent mixing of cloud droplets (i.e., the first right-hand-side term in Eq. (4)) are numerically solved together using an explicit time-stepping method with adaptive step sizes, where Δt in Eqs. (6) and (7) is subdivided into smaller time steps to fulfill the stability requirement.

To aid in the interpretation of some of the results presented later in Sect. 5.2, it is worth clarifying that, when the cloud fraction is non-zero in contiguous layers in the same grid column, the cloud fraction may vary in the vertical direction, and hence cloud bases may be identified in multiple layers according to the criterion stated at the beginning of this section. In this paper, we refer to the cloud base identified in a grid box as the local cloud base of that grid box. If the model layer closest to the Earth's surface is cloud-free, then the interface between the lowest model layer with non-zero cloud fraction and the layer below can be referred to as, in a loose sense, the main cloud base (or lowest cloud base) in the grid column. The DROPMIXNUC parameterization does not distinguish between these two types of cloud base. Both of them lead to the source of cloud droplet number expressed by Eq. (9).

Given the definition of local cloud base noted here, one can imagine that, if a grid column has constant cloud fractions across multiple layers, then these layers will have only a main base and no local base in other cloudy layers. In such a scenario, considering that droplet nucleation by cloud regeneration inside persisting clouds (i.e., the first right-hand-side term in Eq. (5)) is neglected, if the cloud fractions do not increase with time, there will be no droplet nucleation above the main cloud base, and turbulent transport will be the only possible source of the number of cloud droplets in the cloudy layers away from the main base. This design feature is important to keep in mind for some of the discussions in Sect. 5.2.

2.3.4 Droplet evaporation after turbulent transport into clear air

When cloud droplets are transported by turbulence from a layer with a larger cloud fraction into a layer with a smaller cloud, a fraction of the droplets are evaporated, based on the difference in the two cloud fractions, and the cloud-borne aerosols therein are resuspended (Fig. 2d). This is done for transport to both the layer above and the layer below.

As a side note for readers working with EAMv2, we emphasize here that the turbulent mixing of cloud droplet number (and all aerosol species) is numerically solved inside the DROPMIXNUC parameterization and is not handled by CLUBB.

2.4 Implementation of CDNC_min in EAMv2

The lower bound of CDNC, CDNC_min, is set to 10 cm⁻³ in the default EAMv2. Because the model's time integration uses sequential splitting to numerically couple most of the parameterized processes, the value of CDNC in a grid box, if diagnosed, would change multiple times within the 30 min time step depicted in Fig. 1a, as well as during each of the 5 min subcycles used for CLUBB, DROPMIXNUC, and MG2. For our discussions here, it is important to clarify that the lower bound is applied after the DROPMIXNUC parameterization has been calculated and the resulting tendency of cloud droplet number has been applied, and before the cloud microphysical processes are calculated; see boxes III, ^∗, and IV in Fig. 1a.

Due to the aforementioned sequence of calculations, the lower bound, CDNC_min, directly affects all cloud microphysics processes in the MG2 cloud microphysics parameterization that depend on cloud droplet number concentration. In addition, CDNC_min indirectly affects the parameters of the cloud droplet size distribution that are diagnosed near the end of MG2 (box IV in Fig. 1a) and passed subsequently to the radiation parameterization for the calculation of cloud optical properties (box V in Fig. 1a).

3 Numerical experiments and analysis methods

As mentioned in the introduction, this study aims to obtain a process-level characterization of ultra-low CDNCs in stratiform and shallow convective clouds simulated by EAMv2. Our overall strategy is to first identify the geographical locations where ultra-low CDNCs occur most frequently, as well as the locations where the estimated ERF_aer is most sensitive to CDNC_min (Sect. 4). We then examine characteristics of the cloud formation processes associated with ultra-low CDNCs (Sect. 5). Additional sensitivity experiments are used to assist further reasoning about ultra-low CDNCs and their relation to ERF_aer (Sect. 6). Two types of simulation are presented in this paper, as explained next.

Table 1 EAMv2 simulations presented in the paper. “PD” and “PI” refer to the aerosol and precursor emission scenarios, with PD being the 2005–2014 average and PI being values representing the year 1850. The free-running simulations were climatological runs. The nudged simulations were constrained by winds from the ERA-Interim reanalysis. Further details can be found in Sect. 3.

3.1 Nudged simulations for estimating ERF_aer

In this study, ERF_aer was estimated following the methodology used by the Aerosol Comparisons between Observations and Models (AeroCom) intercomparison activities (see, e.g., Gliß et al., 2021; Myhre et al., 2013). Simulations were conducted with interactive atmosphere and land surface models. The sea surface temperature and sea ice extent were prescribed using climatological values. Pairs of otherwise identical simulations were conducted using the present-day (PD) or pre-industrial (PI) emissions of aerosols and precursors, with the PD emissions set to the 2005–2014 averages and the PI emissions set to values representing the year 1850. ERF_aer was calculated as the PD–PI difference in the net radiative flux at the top of the model atmosphere.

These pairs of PD–PI simulations were nudged to the ERA-Interim reanalysis (Dee et al., 2011) to help distinguish signals of aerosol effects from noise caused by natural variability (Kooperman et al., 2012). Horizontal winds in the altitude range of approximately 0.5 to 850 hPa (i.e., model levels 5 to 58 out of a total of 72) were nudged with a relaxation time of 6 h. The reanalysis data were read every 6 h and linearly interpolated to each model time step (Sun et al., 2019; Zhang et al., 2022b). Winds in the near-surface layers, as well as air temperature and humidity, were not nudged, so as to retain good consistency between the climatology of the nudged simulations and the climatology of the free-running simulations (Timmreck and Schulz, 2004; Zhang et al., 2014; Sun et al., 2019; Zhang et al., 2022a, b). Zhang et al. (2022a) showed that 1-year nudged E3SM simulations were sufficient for revealing key signals in the annually averaged global mean, zonal mean, and global patterns of ERF_aer. In this study, however, the sensitivity of ERF_aer to CDNC_min (i.e., ΔERF_aer caused by ΔCDNC_min) was also examined in different seasons (see Fig. 4); this requires more simulation years to distinguish signal from noise. Hence we chose to conduct one pair (PD and PI) of 10-year nudged simulations from 2005 to 2014 using the default EAMv2 (i.e., configuration nc10 in Table 1, with CDNC_min = 10 cm⁻³) and a second pair with CDNC_min set to 0 (i.e., configuration nc00 in Table 1).

Inspired by the results described in Sects. 4 and 5, we conducted several additional pairs of PD–PI simulations. Group II listed in Table 1 had CDNC_min set to 10 cm⁻³ but only in the lower troposphere (between ≈600 hPa and the Earth's surface, configuration nc10_600hPa) or only in grid boxes with a cloud fraction higher than 0.9 (configuration nc10_f0.9). The experiments in group III used no lower bound but either the characteristic updraft velocity used in the droplet nucleation parameterization was artificially enhanced by a factor of 10 (configuration nc00_w10) or both the characteristic updraft and the turbulent mixing coefficient were enhanced by a factor of 10 (nc00_w10k10, Table 1). The group-II and group-III experiments were used for assessing the model's responses in the annual mean ERF_aer; hence, these simulations were integrated for 1 year after 3 months of spin-up.

3.2 Free-running simulation for characterizing ultra-low CDNCs and associated atmospheric conditions

To carry out a process-level analysis of the occurrences of ultra-low CDNCs in EAM, we present a 10-year free-running simulation conducted under forcing conditions of the year 2010 and with the lower bound CDNC_min turned off (i.e., configuration nc00). The reason for not using nudging in this simulation was that free-running simulations were more convenient to set up and would reflect the model's original behavior. That said, we do not expect that nudging horizontal wind will substantially modify the cloud features discussed here.

Our analyses of the free-running simulation relied heavily on the online budget analysis and conditional sampling capabilities developed for EAM by Wan et al. (2022), as well as offline interactive data exploration enabled by the visualization application ParaView. During a simulation, CDNCs were diagnosed at the point in the time loop where ultra-low values would be bounded for the first time in each “main” time step of 30 min when the lower bound was in use, as depicted by the dashed box marked “(*)” in Fig. 1a. CDNCs in the range 0.01–9.99 cm⁻³ are referred to as “ultra-low” in this paper. We note that this range excludes the cases with CDNC=0, as the zeros typically correspond to cloud-free conditions and are not the focus of this study. For comparison with the ultra-low cases, CDNCs of 30–50 cm⁻³ and higher than 50 cm⁻³ were also sampled and are discussed in Sect. 5.

4 Ultra-low CDNCs and ERF_aer

This section aims to assess the frequencies and locations of the occurrences of ultra-low CDNCs in stratiform and shallow convective clouds simulated by EAMv2, as well as to identify the locations at which the simulated ERF_aer is most sensitive to the lower bound, CDNC_min.

Figure 3 (a) Global annual mean ERF_aer simulated with the various configurations of EAMv2 listed in Table 1. (b) Sensitivity of the global annual mean ERF_aer to model configuration, shown as differences with respect to configuration nc00. In both panels, the filled rectangles and the numbers noted therein correspond to the model year of 2011. For the first two rectangles from the top in panel (a) and the first rectangle from the top in panel (b), the whiskers visualize the ranges of 1 year annual mean among the 10 simulation years, and the three numbers noted in each pair of parentheses are the minimum, mean, and maximum values of the 1-year annual averages.

Figure 4 (a) 10-year mean annually averaged ERF_aer, diagnosed from a pair of PD and PI simulations nudged to the meteorology of the years 2005–2014, conducted using EAMv2 with no lower bound for CDNC (i.e., configuration nc00). (b) Changes in 10-year annual mean ERF_aer caused by imposing a lower bound of CDNC_min=10 cm⁻³ (configuration nc10 minus nc00). (c) As panel (b) but showing the relative changes. In panels (a–c), masked out in white are the geographical locations where the ERF_aer in configuration nc00 is statistically small (see Sect. 4.1). (d–f) 10-year mean zonal mean monthly ERF_aer in the Northern Hemisphere: (d) in configuration nc00, (e) in configuration nc10, and (f) the differences. Similar figures showing the shortwave (SW) component of ERF_aer and the SW component of the aerosol indirect effect (AIE) can be found in Figs. S2 and S3. Further details and discussions can be found in Sect. 4.1.

4.1 Sensitivity of ERF_aer to CDNC_min

We first compare the ERF_aer simulated by EAMv2 with or without using the lower bound. The 10-year mean global and annual mean ERF_aer are −1.71 W m⁻² in the nc00 configuration and −1.43 W m⁻² in the nc10 configuration (Fig. 3a). Figure 4a shows a global map of the annual mean ERF_aer derived from the nc00 configuration (i.e., no lower bound). Figure 4b and c show the signed absolute and relative changes, respectively, caused by imposing a CDNC_min of 10 cm⁻³. In all three panels, masked out in white are the geographical locations with statistically small ERF_aer in nc00, i.e., where an absolute value of the 10-year annual mean ERF_aer is smaller than the estimated 2 standard deviations of the 10-year annual mean.

A comparison between Fig. 4a and b suggests that the lower bound reduces ERF_aer globally except over the major deserts. In the middle latitudes of both hemispheres, the spatial pattern of the change (Fig. 4b) largely matches the spatial pattern of ERF_aer itself (Fig. 4a), meaning that there is a rough proportionality between the magnitudes of change and the baseline values.

In the Arctic, however, the annual mean absolute changes are large, often exceeding 1 W m⁻² (Fig. 4b), while the baseline values are moderate (within −2 W m⁻² in general, Fig. 4a). Therefore, this region features the strongest relative reduction in the magnitude of ERF_aer, of the order of 50 % to 100 %, when a CDNC_min of 10 cm⁻³ is imposed (Fig. 4c). This strong sensitivity in the annual mean ERF_aer is dominated by the changes in boreal summer, especially June and July, which can be seen in the zonally averaged monthly mean ERF_aer differences shown in Fig. 4f. In the Northern Hemisphere (NH) middle latitudes, the 10-year mean monthly zonal mean changes are typically within 1.5 W m⁻², while the changes near the North Pole are 5 W m⁻² or larger in June and July.

As can be expected from the results reported in Zhang et al. (2022a), the changes in ERF_aer seen in Fig. 4 are dominated by the changes in the shortwave ERF_aer (see Fig. S2), which, in turn, are attributable primarily to the changes in the shortwave component of the aerosol indirect effect (see Fig. S3).

Figure 5 (a–c) 10-year mean annually averaged fraction of cloudy time steps and model layers that exhibit ultra-low CDNCs, i.e., the conditional probability $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\mathrm{|}\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$. (d–f) 10-year mean annually averaged fraction of all model time steps and model layers that are considered “cloudy” because the grid box mean cloud liquid mass mixing ratio of stratiform and shallow convective clouds, $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$, is higher than a threshold q_l,min, i.e., the probability $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$. The three columns correspond to three different choices of q_l,min used in conditional sampling: q_l,min = 10⁻¹⁵ g kg⁻¹ in the left column, q_l,min = 10⁻⁵ g kg⁻¹ in the middle column, and q_l,min = 0.02 g kg⁻¹ in the right column. The fractions (probabilities) were diagnosed in the 10-year free-running nc00 simulation explained in Table 1, i.e., EAMv2 with no lower bound for CDNC. The geographical distribution of the joint probability $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\cap \stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ can be found in Fig. S4. Further details and discussions can be found in Sect. 4.2.

4.2 Geographical distribution of ultra-low CDNCs

Figure 5a is a global map of the 10-year mean annually averaged fraction of cloudy time steps and model layers that exhibit ultra-low CDNCs in the nc00 simulation. Here, “cloudy” is defined as $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}}$ = 10⁻¹⁵ g kg⁻¹, as is the case at the point in the model's time loop where CDNC_min is applied, and $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$ is the instantaneous grid box mean cloud liquid mass mixing ratio in stratiform and shallow convective clouds. All vertical layers (72 in total in EAMv2) were included in the calculation of the fraction. Over the mid- and high-latitude oceans in both hemispheres, 25 % to 55 % of the cloudy time steps and vertical levels are found to have ultra-low CDNCs, and we see trends of higher percentages at higher latitudes.

The fraction shown in Fig. 5a can be interpreted as a conditional probability, $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\mathrm{|}\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ and was diagnosed using

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\mathrm{|}\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})\\ & =P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\cap \stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})/P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}}).\end{array}\end{array}$$

To avoid drawing undue attention to large values of the conditional probability resulting from singularities associated with very small values of the denominator $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, Fig. 5d shows $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, namely the average fraction of cloudy cases among all model time steps and all vertical levels. Compared with the low latitudes and most continental areas, the mid- and high-latitude oceans are associated with a more frequent presence of cloud liquid mass. Hence, the high frequencies of ultra-low CDNC depicted by Fig. 5a over the mid- and high-latitude oceans are meaningful signals. The geographical distribution of the joint probability $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\cap \stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ can be found in Fig. S4a.

One might also wonder whether the high frequencies seen in Fig. 5a and d have resulted from using a very low threshold of q_l,min = 10⁻¹⁵ g kg⁻¹ to define “cloudy” (i.e., cloud liquid present). Since $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$ is the grid box mean value and the grid boxes are large (approximately 200 km × 200 km), and since numerical discretization and floating-point calculations both introduce error, we can imagine a simulation producing $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$ values that exceed 10⁻¹⁵ g kg⁻¹ but are nevertheless very small from a practical perspective and hence unlikely to have large impacts on the mean climate. To avoid drawing undue attention to such cases, we diagnosed, in the same nc00 simulation, fractional frequencies similar to those shown in Fig. 5a and d but with q_l,min set to 10⁻⁵ g kg⁻¹ in conditional sampling. The resulting global maps are presented in Fig. 5b and e, which turn out to closely assemble the results obtained using the smaller threshold of 10⁻¹⁵ g kg⁻¹. Further increasing q_l,min to 0.02 g kg⁻¹ leads to substantial decreases in the count of cloudy cases (Fig. 5f versus Fig. 5e), but the mid- and high-latitude oceans still stand out as the regions with most-frequent cloud liquid presence (Fig. 5f). Overall, the three sets of diagnostics presented in the figure suggest that, over the middle- and high-latitude oceans, 25 % to 55 % of the cases (i.e., grid boxes and time steps) with cloud liquid presence are associated with ultra-low CDNCs. The feature is robust across a wide range of q_l,min used for defining cloud liquid presence.

Between the two hemispheres, the Southern Hemisphere (SH) shows slightly more frequent occurrences of ultra-low CDNCs (Figs. 5a–c and S4a–c). The absolute changes in ERF_aer caused by imposing CDNC_min, however, are substantially smaller in the SH than in the NH (Fig. 4b). This is not surprising, given the significantly lower aerosol concentrations and their PD–PI differences in the SH. Since this study was originally motivated by the ERF_aer in EAMv2, the analyses in the remainder of the paper will focus on the NH.

In southeast (SE) China and the SE United States, about 15 % to 25 % of the time steps and vertical levels with cloud liquid presence are associated with ultra-low CDNCs, regardless of the value of $q_{\mathrm{l},min}$ used in conditional sampling (Fig. 5a–c). These percentages might seem unexpectedly high for continental regions with generally high aerosol concentrations. In the following, these two continental regions are further analyzed, together with the Arctic region and the North Pacific storm track. The latter two regions were chosen because they feature both high frequencies of occurrence of ultra-low CDNCs (Fig. 5a–c) and strong sensitivities of ERF_aer to CDNC_min (Fig. 4b).

Figure 6 Regionally averaged vertical profiles of the annual mean fractional frequency of occurrence of various conditions. As in Fig. 5, the upper row (a–d) shows the conditional probability $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\mathrm{|}\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, i.e., the fraction of cloudy cases that exhibit ultra-low CDNCs, and the lower row (e–h) shows the probability $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, i.e., the fraction of model time steps and vertical layers that are considered “cloudy” because $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}}$. The four columns correspond to different regions; from left to right are the North Pacific storm track, the Arctic region, SE China, and SE United States, respectively. Lines are 10-year averages, and color shading shows 2 standard deviations of ten 1-year averages. The solid, dashed, and dotted lines correspond to different q_l,min values used in conditional sampling. The fractions (probabilities) were diagnosed from the 10-year free-running nc00 simulation explained in Table 1 (i.e., EAMv2 with no lower bound for CDNC). Vertical profiles of the joint probability $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\cap \stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ can be found in Fig. S5. Further details and discussions can be found in Sect. 4.3.

4.3 Vertical distribution of ultra-low CDNCs

Based on the findings in Sects. 4.1 and 4.2, we selected four focus regions, namely the North Pacific storm track (oceanic areas in 30–60° N, 120–240° E), the Arctic region (all geographical locations northward of 66.5° N), SE China (land areas in 25–35° N, 110–120° E), and the SE United States (land areas in 30–40° N, 270–280° E), to examine the vertical distribution of the ultra-low CDNCs. Similar to Fig. 5, Fig. 6a–d show the 10-year mean fractional frequency of occurrences of ultra-low CDNCs in samples with cloud liquid presence, i.e., $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\mathrm{|}\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, and Fig. 6e–h show the fraction of samples featuring $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}}$ among all time steps and grid boxes in the region and vertical layer, i.e., $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$. The vertical profiles of the joint probability $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\cap \stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ can be found in Fig. S5. The four columns in Fig. 6 show results for the four regions. The solid, dashed, and dotted lines in each panel are 10-year mean annual averages corresponding to q_l,min = 10⁻¹⁵, 10⁻⁵, and 0.02 g kg⁻¹, respectively. Two standard deviations of 1-year averages are shown as color shading to indicate interannual variability.

The four regions differ in their frequencies of cloud liquid presence, as expected (Fig. 6e–h), but similar features can be seen in the fractions of cloudy cases exhibiting ultra-low CDNCs (Fig. 6a–d). In the upper troposphere (with pressure values lower than 300 or 400 hPa, depending on the region), the occurrences of ultra-low CDNCs in liquid-present cases vary substantially with q_l,min, possibly due to the dominance of ice clouds and very small values of $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ at those altitudes. In the middle and lower troposphere (pressure>400 hPa), the fractions of liquid-present cases associated with ultra-low CDNCs are rather insensitive to q_l,min (Fig. 6a–d). Compared with the 10-year averages, the interannual variabilities are small, suggesting that the mean frequencies shown here are statistically significant.

The clouds in EAMv2 in the upper troposphere are often cirrus, while those occurring between 600 hPa and the surface are typically mixed-phase or warm (liquid) clouds. To find out which of these cloud types play more important roles in causing the sensitivity of ERF_aer to CDNC_min, we conducted the sensitivity experiment nc10_600 hPa listed in Table 1, in which CDNC_min = 10 cm⁻³ was applied only to the model's vertical layers 48 to 72, corresponding to pressure values equal to or higher than about 600 hPa when the surface pressure was 1000 hPa. Compared with nc00 (no lower bound), bounding CDNC in all layers (i.e., nc10) reduces the magnitude of the global annual mean ERF_aer by 0.34 W m⁻² in the year 2011 and by 0.19–0.36 W m⁻² in the 10-year period of 2005 to 2014 (Fig. 3b). Bounding CDNCs between 600 hPa and the surface gives a reduction of 0.33 W m⁻² in the year 2011, suggesting that the ultra-low CDNCs in the lower troposphere play a dominant role in determining the sensitivity of ERF_aer to CDNC_min.

For the lower-tropospheric layers between 600 and about 900 hPa, Fig. 6b reveals that about 40 % of the cloudy cases (i.e., cloudy time steps and grid boxes) in the Arctic region are associated with ultra-low CDNCs. Over the North Pacific storm track (Fig. 6a), SE China (Fig. 6c), and the SE USA (Fig. 6d), the typical percentages are between 20 % and 40 %. These are rather high annual mean percentages, especially for the polluted regions.

4.4 Cloud morphology

Having identified the locations of the most frequent occurrences of ultra-low CDNCs, we now start to examine the features of the model atmosphere in those cases, focusing on the lower troposphere and on the four regions discussed earlier (Fig. 6).

Figure 7 (a–d) Bivariate histograms of the grid box mean cloud liquid mass mixing ratio, $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$, and cloud fraction, f_c, for the ultra-low CDNCs occurring in model layers with nominal pressure higher than 600 hPa. (e–h) Bivariate histograms of the in-cloud droplet mass concentration (CDMC) and f_c. The model output used in this analysis was 120 d of 6-hourly instantaneous values from January, April, July, and October (30 d per month) of the first simulation year after spin-up. As in Fig. 6, results from different focus regions are shown in different columns, as indicated above each plot panel. Further details and discussions can be found in Sect. 4.4.

Our preliminary examination of the grid box mean cloud liquid mass mixing ratio ($\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$, unit: g kg⁻¹), the in-cloud droplet mass concentration (CDMC, unit: g cm⁻³), and cloud fraction (f_c, unitless) suggested that the cases with ultra-low CDNCs showed moderate to high composite mean values for all three quantities. A subsequent investigation using instantaneous model output revealed the existence of two distinct regimes. Figure 7a–d are bivariate histograms of instantaneous $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$ and f_c in the four focus regions, presented for the ultra-low CDNC cases sampled in the lower troposphere from 120 d of 6-hourly instantaneous model output. Figure 7e–h are bivariate histograms of CDMC and f_c. To include different seasons, the 120 d used in this analysis consisted of 30 consecutive days per month from each of January, April, July, and October of the first simulation year. The results from the four geographical regions turn out to show very similar features: the majority of the samples with ultra-low CDNCs are associated with f_c>0.9, typical values of $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$ between 0.01 and 1 g kg⁻¹, and in-cloud droplet mass concentrations from 0.01 to 1 g cm⁻³. There is a second regime featuring very small cloud fractions (f_c<0.1), typical values of $\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}$ lower than 0.01 g kg⁻¹, and in-cloud droplet mass concentrations spanning a wide range, from about 10⁻⁵ to 1 g cm⁻³. Only a very small portion of the samples is associated with cloud fractions between 0.1 and 0.9.

These results naturally trigger the following questions:

1. Which of the two regimes plays a larger role in affecting the simulated global mean ERF_aer?

2. Where do the samples in each regime occur in terms of cloud morphology? In other words, do the ultra-low CDNCs associated with small or large cloud fractions typically occur inside a cloudy region or at the boundaries between cloudy and cloud-free regions? Do they occur in small and scattered cloudy regions or are they associated with large (e.g., synoptic-scale) cloud systems?

Figure 8 Vertical cross-sections of (a, d) cloud fraction, (b, e) in-cloud CDNC, and (c, f) in-cloud droplet mass concentration along the prime meridian northward of 66.5° N at 00:00 UTC on 3 July of the first year in the free-running nc00 simulation. The color shading is the same in both rows, but the overlaid quadrilateral frames in the upper row indicate model grid boxes associated with ultra-low CDNCs and small cloud fractions (<0.1), while the quadrilateral frames in the lower row indicate model grid boxes associated with ultra-low CDNCs and large cloud fractions (>0.9). An animation showing vertical cross-sections in the 30 d period from 28 June to 27 July can be found in the Video Supplement. See Sect. 4.4 for further details and discussions.

The answers to these questions may provide important clues for further investigations. For example, when ultra-low CDNCs are associated with very small cloud fractions and are sporadically located (e.g., at the edges of cloudy regions), one can suspect they might have been caused by numerical artifacts of spatial discretization. In contrast, if impactful ultra-low CDNCs are found in large contiguous areas and at specific types of location (e.g., in the center or at the bottom of a synoptic-scale cloudy region), then the ultra-low values may have resulted from certain inherent features of the physical assumptions in the model.

Our sensitivity experiment nc10_f0.9 was designed to answer the first question. In this pair of PD and PI simulations, the lower bound CDNC_min was applied only to grid boxes with cloud fractions higher than 0.9. The change in the year-2011 mean global mean ERF_aer relative to nc00 amounts to a reduction of 0.22 W m⁻² in magnitude, explaining a major portion (about 65 %) of the reduction of 0.34 W m⁻² in magnitude obtained by eliminating all ultra-low CDNCs (Fig. 3b). We also conducted an additional experiment, bounding CDNCs only in grid boxes with cloud fractions lower than 0.1; the resulting global mean ERF_aer turned out to be very similar to that of the nc00 configuration. It then follows that the aggregate effect of the cases with cloud fractions between 0.1 and 0.9 is non-negligible, despite the relatively low case counts in those cloud fraction bins. How the elimination of ultra-low CDNCs affects ERF_aer and its various components in the high-cloud-fraction and medium-cloud-fraction regimes is an interesting topic for future investigations.

Figure 9 (a) Vertically integrated low-level cloud fraction at 00:00 UTC on 8 March of the first model year. The dashed square is the SE China region that panels (b) and (c) zoom into. (b) The semi-transparent color shading shows the cloud fraction values at the point of the model's time loop where the CDNCs are diagnosed and the lower bound is applied. The solid 3D boxes with magenta outlines are the grid boxes featuring ultra-low CDNCs and large cloud fractions (>0.9). (c) As (b), but the solid 3D grid boxes are associated with ultra-low CDNCs and small cloud fractions (<0.1). An animation showing 3D visualizations in the 30 d period from 28 February to 29 March can be found in the Video Supplement. See Sect. 4.4 for further details and discussions.

The second question asked here is not straightforward to answer using numerical metrics, so we chose to resort to visual inspection using ParaView. Given the finding, presented in Sect. 4.1, that ERF_aer simulated by EAMv2 features the strongest sensitivity to CDNC_min in Arctic summer, Fig. 8 presents various vertical cross-sections in the Arctic circle along the prime meridian, captured at 00:00 UTC on 3 July of the first simulation year. The two rows in the figure present the same quantities in the format of color shading: cloud fraction in Fig. 8a and d, CDNC in Fig. 8b and e, and CDMC in Fig. 8c and f. The quadrilateral frames in Fig. 8a–c indicate grid boxes with ultra-low CDNCs and small cloud fractions (f_c<0.1). These cases are found sporadically (Fig. 8a) and are associated with low in-cloud droplet mass concentrations (Fig. 8c). This combination of features seems physically plausible. The quadrilateral frames in Fig. 8d–f indicate EAMv2's grid boxes with ultra-low CDNCs and large cloud fractions (f_c>0.9). These grid boxes are located in the main body of a multi-layer cloud system (Fig. 8d), where CDNCs are systematically low across a number of vertical layers and grid columns (Fig. 8e). The in-cloud droplet mass concentrations in these grid boxes, however, can be moderate or high (e.g., more than 0.05 g m⁻³, Fig. 8f), which is consistent with the bivariate histogram shown in Fig. 7f. This combination of features, i.e., ultra-low CDNC but high cloud fraction and substantial in-cloud droplet mass concentration, is counterintuitive. While Fig. 8 shows only one time instance, the distinct features of the two types of ultra-low CDNC case are frequently seen in the 30 d of instantaneous output that we examined from June 28 to July 27 of the first simulation year. An animation showing vertical cross-sections for the 30 d period can be found in the Video Supplement.

A similar analysis of cloud morphology has been carried out for the SE China region. To ease the search for relevant cases, we examined the monthly mean regionally averaged vertical profiles of $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\mathrm{|}\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, $P(\text{CDNC\hspace{0.17em}\hspace{0.17em}ultra-low}\cap \stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$, and $P(\stackrel{\mathrm{‾}}{q_{\mathrm{l}}}⩾q_{\mathrm{l},\text{min}})$ in the first simulation year and selected March as the focus month, as all three fractions were relatively high in that month (see third column of Fig. S6). We then selected a model output time, namely 00 UTC on March 8, when the region of 25–35° N, 110–120° E was fully covered by low-level clouds, as judged from the vertically integrated low-level cloud fraction (Fig. 9a). In Fig. 9b and c, cloud fractions in the region of 25–35° N, 110–120° E are shown in semi-transparent color shading, and the grid boxes with ultra-low CDNCs and large (>0.9, Fig. 9b) or small (<0.1, Fig. 9c) cloud fractions are shown as light gray or dark blue solid 3D boxes, respectively, with the colors indicating the cloud fraction values. The two panels, Fig. 9b and c, reveal a feature similar to that seen for the Arctic: ultra-low CDNCs associated with large cloud fractions are found in contiguous grid boxes inside large cloudy regions, while the ultra-low CDNCs associated with small cloud fractions are scattered in space. An animation showing the same kind of 3D visualization for the 30 d period from 28 February to 29 March of the first simulation year can be found in the Video Supplement, and suggests that the distinct features of the two types of ultra-low CDNC case are typical over SE China.

Since the ultra-low CDNCs associated with large cloud fractions have a major impact on the global mean ERF_aer and are intriguingly located inside large cloudy regions, we focus on this category of cases in the next section, using composite analysis to understand characteristics of the droplet formation and removal processes.

5 Cloud droplet budget analyses

Generally speaking, ultra-low CDNCs in stratiform and shallow convective clouds could be caused by weak sources or strong sinks of the number of droplets, or both. To help identify the culprits, we present in Sect. 5.1 a budget analysis for the main groups of droplet formation, transport, and removal processes considered in EAMv2, namely the dynamical core and parameterizations depicted by gray boxes in Fig. 1a and described in Sect. 2.1. After that, Sect. 5.2 zooms into the DROPMIXNUC parameterization to analyze the processes depicted in Fig. 2 and described in Sect. 2.3. All results presented in this section are based on the 10-year free-running simulation using the nc00 configuration, i.e., EAMv2 with no lower bound for CDNC.

Figure 10 Results from the 10-year free-running simulation using the nc00 configuration, i.e., EAMv2 with no lower bound for CDNC. (a) Regionally averaged Arctic July mean frequency of occurrence of grid boxes and model time steps with cloud fraction larger than 0.9 and CDNC (unit: cm⁻³) in three different ranges. (b–d) Composite mean of the grid box mean cloud droplet number tendencies caused by resolved transport, DROPMIXNUC, and cloud microphysics, respectively. (e) Composite mean cloud fraction. (f–h) Composite mean of the grid box mean cloud droplet mass tendencies caused by resolved transport, CLUBB, and cloud microphysics, respectively. Different marks and colors correspond to different CDNC ranges. The marks and lines show the averages of the month of July of 10 simulation years. Color shading indicates 2 standard deviations of the July average of individual years. Further details and discussions can be found at the beginning of Sect. 5 and in Sect. 5.1.

To put the numbers associated with ultra-low CDNCs into context, two additional CDNC ranges, 30–50 cm⁻³ and higher than 50 cm⁻³, were sampled for comparison, all under the condition of cloud fraction being higher than 0.9. For clarification, we note that the conditional sampling was done pointwise (i.e., grid-box-wise rather than column-wise), as we intended to examine the local conditions associated with different CDNC ranges. At any time step, if one grid column had, for example, 2 grid boxes with ultra-low CDNCs, and another grid column had 10 grid boxes with ultra-low CDNCs, these grid boxes would each contribute to the count of ultra-low CDNCs in their respective model layers, and the atmospheric conditions (or process rates) in those grid boxes were recorded to diagnose composite mean values.

The fractional frequency of occurrence of the three CDNC ranges with respect to all model time steps in the month of the year and all grid boxes in the model layer and region is shown in Fig. 10a for Arctic July as a function of altitude (pressure). The marks and lines indicate 10-year averages of the focus month. Color shading indicates 2 standard deviations of the July averages of 10 different years. Figure 10a shows that, in Arctic July, ultra-low CDNCs occur significantly more often than higher concentrations, with a clear trend of higher frequencies of ultra-low concentrations toward the surface. In the figures presented in the remainder of this section, composite mean vertical profiles of various physical quantities are shown. (The results are masked out, i.e., not shown for the specific altitude or sampling condition, if the fractional frequency of occurrence of the condition is lower than 0.1 % in the focus month in any of the 10 years, to avoid showing results for very small samples.)

The composite mean cloud fractions of the sampled CDNC ranges are shown in Fig. 10e for Arctic July. The mean values are generally higher than 0.95 for all three CDNC ranges and are about 0.99 or higher for the ultra-low CDNCs. Due to these very high cloud fractions, although all process rates (and other physical quantities) shown in the remainder of this section are grid box mean values, we expect the in-cloud mean values to be similar.

5.1 Source problem or sink problem?

Figure 10b–d present the composite mean cloud droplet number tendencies in Arctic July caused by the resolved transport, DROPMIXNUC, and cloud microphysics. The parameterized vertical transport by updrafts and downdrafts in deep convection and the detrainment from deep convection have negligible magnitudes and hence are not shown. This budget analysis shows that, overall, DROPMIXNUC is the primary source of cloud droplet number in this region and month of the year, while cloud microphysics and resolved transport are both significant sinks.

A comparison among the three CDNC ranges reveals that both the sources and the sinks associated with ultra-low CDNCs are markedly weaker (closer to zero) than those seen for higher CDNCs (Fig. 10b–d). This feature is also seen in the other three regions (see Figs. S7–S9), suggesting that the ultra-low CDNCs are probably caused by weak sources, rather than strong sinks.

5.2 Droplet nucleation and turbulence

Figure 10f–h present a budget analysis similar to Fig. 10b-d but for the mass of cloud droplets. As expected from the construct of EAMv2 and its parameterizations, in most model layers in the lower troposphere (except for the very few layers closest to the Earth's surface), the turbulence and cloud parameterization CLUBB acts as the primary source of droplet mass, while cloud microphysics is the primary sink. Somewhat counterintuitively, when comparing the three CDNC ranges, we see that the cases with ultra-low CDNCs are associated with substantially stronger mass sources and sinks (Fig. 10f to h), despite their much weaker number sources and sinks revealed by Fig. 10b–d. This contrast in tendencies is consistent with the feature revealed earlier in Fig. 8, that the grid boxes with ultra-low CDNCs and large cloud fractions are often associated with medium to high in-cloud droplet mass concentrations. Furthermore, the contrast between strong mass sources/sinks and weak number sources/sinks is also seen clearly in the other three regions, see Figs. S7–S9.

Figure 11 Regionally averaged composite mean of various quantities in Arctic July for three different CDNC ranges (CDNC unit: cm⁻³). (a) Grid box mean droplet number mixing ratio tendency due to aerosol activation at the local cloud base in the grid box (Eq. (9) and Sect. 2.3.3). (b) Grid box mean droplet number mixing ratio tendency due to aerosol activation associated with cloud fraction increase (Eq. (7) and Sect. 2.3.2). (c) Grid box mean droplet number mixing ratio tendency due to turbulent mixing. (d) Fractional area of the local cloud base in the grid box, i.e., f_b in Eq. (9), defined at the beginning of Sect. 2.3.3. (e) Characteristic updraft velocity used for calculating aerosol activation, see Eq. (2). (f) Subgrid variance of vertical velocity, see Eq. (3). (k) Skewness of the subgrid vertical velocity calculated by CLUBB. All results were derived from the 10-year free-running simulation using configuration nc00, i.e., EAMv2 with no lower bound for CDNC. The marks and lines show the averages of the month of July of 10 simulation years. Color shading indicates 2 standard deviations of the July average of individual years. Further details and discussions can be found in Sect. 5.2.

To investigate why the droplet number sources associated with ultra-low CDNCs are weak while the corresponding droplet mass sources are strong, we decomposed the number tendency caused by the DROPMIXNUC parameterization, i.e., the quantity shown in Fig. 10c, into the various contributors described in Sect. 2.3: droplet nucleation at the local cloud base in the grid box (Fig. 11a), droplet nucleation associated with cloud fraction increase in the grid box (Fig. 11b), and turbulent mixing (Fig. 11c). The evaporation associated with cloud fraction decrease and the evaporation of droplets detrained by turbulence are both negligible in magnitude and hence not shown. The figure shows that the cases with higher CDNCs (30–50 cm⁻³ or >50 cm⁻³) are associated with significant droplet number sources due to nucleation at the local cloud base in the grid box (Fig. 11a); positive tendencies are seen in almost all layers in the lower troposphere and peaks of nucleation rate are found in the near-surface layers (Fig. 11a, green and blue). The tendencies associated with turbulence mixing are negative (Fig. 11c, green and blue), suggesting that droplets are transported away from where they are nucleated. Nucleation associated with cloud fraction increase is negligible in most layers with pressure higher than about 800 hPa. The positive net droplet number tendencies seen in Fig. 10c for the DROPMIXNUC parameterization in these higher CDNC ranges turn out to be small residuals between cloud-base nucleation (a substantial source) and turbulent mixing (a substantial sink). In contrast, in the cases with ultra-low CDNCs, the droplet number tendencies are very close to zero for all three processes (Fig. 11a–c, red) and for the total tendency caused by DROPMIXNUC (Fig. 10c, red).

Figure 11a reveals near-surface peaks of cloud-base droplet nucleation in the cases of higher CDNCs but negligible nucleation associated with ultra-low CDNCs. In Sect. 2.3.3, we explain that this source of cloud droplet number is associated with the turbulent influx of activated aerosols at the local cloud base (i.e., the cloud base in the current grid box) and that the resulting droplet number tendency is proportional to the fractional area of the local cloud base (i.e., f_b in Eq. (9)). Figure 11d shows that, when ultra-low CDNCs are simulated (red line), the fractional area of the local cloud base, f_b, has composite mean values very close to 0, which explains, at least partially, the lack of droplet nucleation in such cases: cloud-base nucleation cannot occur when there is no cloud base.

Section 2.3.3 also points out that, when there is neither droplet nucleation at the local cloud base nor nucleation associated with cloud fraction increase, there still can be, in principle, a source of droplet number caused by turbulent mixing. This potential source, however, also appears to be absent in the cases with ultra-low CDNCs (Fig. 11c, red line). We note that a significant source of droplet number caused by turbulent mixing would require (I) substantial nucleation in other grid layers to cause vertical gradients and (II) sufficiently strong turbulent eddies to carry out down-gradient transport. The online conditional sampling performed in this study was pointwise (grid-box-wise), as the current software infrastructure cannot easily capture column-wide information for each of the sampled grid boxes. Therefore we do not have good model output to reason about factor (I). As for factor (II), Fig. 11e shows the composite mean of $w_{\text{act}}^{\ast }$ defined in Eq. (2) and Fig. 11f shows the σ_w defined in Eq. (3). While $w_{\text{act}}^{\ast }$ is used for calculating the activated aerosol fractions and has an assumed minimum of 0.1 m s⁻¹, as in various other models (see, e.g., discussions in Poku et al., 2021), σ_w does not have this assumed minimum and is an indicator of the strength of turbulence. In the cases of ultra-low CDNCs in Arctic July, the mean σ_w is about 0.2 m s⁻¹ in the lowest model layer but quickly drops to 0.04 m s⁻¹ or less at about 950 hPa and higher altitudes (Fig. 11f), and the values are systematically and substantially smaller than those seen for the higher CDNC ranges (Fig. 11f). This means that, even though droplet nucleation might happen in other layers of the grid column, turbulent transport from those layers is expected to be weak as turbulence is weak. Overall, it is understandable that when the simulated cloudy regions are relatively deep (i.e., cover multiple grid layers) and the cloud fractions vary little in the vertical direction as well as in time, the lack of local nucleation and the rather weak turbulence mixing can lead to a lack of droplet number sources in model layers away from the main cloud base.

Counterparts of Fig. 11 showing the results for January over the North Pacific storm track, March over SE China, and January over the SE United States can be found in Figs. S10–S12. The focus months were chosen for their relatively high frequencies of occurrence of ultra-low CDNCs (see Fig. S6). All the examined regions show the same qualitative features: compared with the cases with higher CDNCs, the ultra-low CDNCs are associated with negligible droplet nucleation, lack of local cloud base, weak turbulence, and negligible turbulent mixing. These results, especially the fact that all sources and sinks represented by the DROPMIXNUC parameterization are negligibly small when ultra-low CDNCs are found, provide additional evidence to support the conclusion drawn at the end of Sect. 5.1 that the ultra-low CDNCs in EAMv2 are probably caused by weak sources rather than strong sinks.

Figure 12 Regionally averaged monthly mean composite mean of cloud condensation nuclei (CCN) number concentration at 0.1 % supersaturation in different regions and focus months: (a) January over the North Pacific storm track, (b) Arctic July, (c) March over SE China, (d) January over the SE United States. The marks and lines show the averages of the focus month over 10 simulation years. Color shading indicates 2 standard deviations of the monthly average of different years. Further details can be found in Sect. 5.3.

5.3 Droplet nucleation problem or aerosol concentration problem?

At this point, one might start wondering whether the lack of cloud droplet number sources identified in the cases of ultra-low CDNCs is caused by the lack of aerosol particles (APs), noting that, if APs were abundant, even small activation fractions and relatively weak turbulent transport might be sufficient to avoid ultra-low CDNCs. To help answer this question, we present in Fig. 12 the composite mean number concentrations of cloud condensate nuclei (CCN) at 0.1 % supersaturation. The CCN number concentrations (unit: cm⁻³) were diagnosed during the simulation using

$$\begin{array}{}\text{(11)}& {\mathit{\rho }}_{\text{air}}{\stackrel{\mathrm{‾}}{N}}_{\text{CCN}}={\mathit{\rho }}_{\text{air}}\sum _{i=\mathrm{1}}^M\left[f_{\mathrm{a},i}^{\ast }\left({\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}+{\stackrel{\mathrm{‾}}{N}}_{\mathrm{c},i}\right)\right],\end{array}$$

where ρ_air is air density and ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}$ and ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{c},i}$ are the grid box mean interstitial and cloud-borne aerosol number mixing ratios, respectively. The activation fractions, $f_{\mathrm{a},i}^{\ast }$, $i=\mathrm{1},\mathrm{\dots },M$, were calculated following Eqs. (13) and (15) in ARG2000, namely,

$$\begin{array}{}\text{(12)}& f_{\mathrm{a},i}^{\ast }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}-\text{erf}\left(u_i^{\ast }\right)\right],\end{array}$$

with

$$\begin{array}{}\text{(13)}& u_i^{\ast }={\displaystyle \frac{\mathrm{2}\ln (S_{\mathrm{m},i}/S^{\ast })}{\mathrm{3}\sqrt{\mathrm{2}}\ln {\mathit{\sigma }}_i}}.\end{array}$$

Here S_m,i is the critical supersaturation of APs with the geometric mean radius of aerosol mode i (Eq. (9) in ARG2000 and Eq. (8) in Abdul-Razzak et al., 1998), σ_i is the geometric standard deviation of aerosol mode i (Eq. (2) in ARG2000), and S^∗ is the prescribed supersaturation of 0.1 % instead of the maximum supersaturation, S_max, in ARG2000. In other words, the diagnosed CCN concentration is the number concentration of APs of all sizes and compositions in the grid box that would be activated at 0.1 % supersaturation.

In Fig. 12, the composite mean CCN concentrations are shown for the three CDNC ranges and for all four regions discussed. The thin black vertical line in each panel is a reference line corresponding to a number concentration of 10 cm⁻³. For all regions and altitudes shown in the figure, the composite mean CCN concentrations associated with ultra-low CDNCs are higher than 10 cm⁻³. This is especially so in SE China, where the mean CCN concentrations between 700 hPa and the surface are generally higher than 400 cm⁻³. In other words, if all APs in the grid box were deemed available for growth into cloud droplets, and if the balance between the supersaturation production in the ambient atmosphere and the supersaturation consumption by the APs' size growth translated to an S_max of 0.1 % in the ARG2000 aerosol activation parameterization, we would expect cloud droplet number concentrations about 2 orders of magnitude higher than the simulated values. It follows that the calculation of droplet nucleation (rather than the simulated AP concentration) plays a primary role in causing ultra-low CDNCs in this region and time of year.

The composite mean CCN concentrations associated with ultra-low CDNCs in the lower troposphere turn out to be in the range of about 20–50 cm⁻³ in January over the North Pacific storm track, 30–40 cm⁻³ in Arctic July, and 10–120 cm⁻³ in January over the SE United States. These concentrations are higher than 10 cm⁻³ but not as dramatically as in SE China. Also, noting that in these regions the CCN concentrations associated with different CDNC ranges are significantly different (i.e., the 2-standard-deviation ranges hardly overlap), we speculate that the AP concentrations also play a role in causing ultra-low CDNCs in these regions. Future work evaluating both the CDNC and AP concentrations against observational data will be helpful. It will be useful to carry out a follow-up study to understand the reasons for the differences between the diagnosed CCN concentrations and the simulated CDNC concentrations, identifying which of these differences in various regions and in different CDNC ranges are physically justifiable and which ones are indicative of the need for model improvement.

6 Additional experiments and discussions

The results in the previous section suggest that weak droplet nucleation is an important cause of ultra-low CDNCs in stratiform and shallow convective clouds in EAMv2. We now present some additional investigations and discussions along this line.

6.1 Motivation and design

Since the results presented in preceding sections suggest that the sampled cases with ultra-low CDNCs and large cloud fractions feature weak turbulence, one might ask the following questions:

• How would the simulated CDNCs and global mean ERF_aer respond if stronger subgrid updraft velocities were simulated by EAMv2 to cause stronger aerosol activation?

• How would the simulated CDNCs and global mean ERF_aer respond if stronger turbulence were simulated by EAMv2 to cause not only more aerosol activation but also stronger subgrid transport of droplets from the cloud base? In other words, would it be worth adopting the assumption of vertically constant CDNC across contiguous cloudy grid boxes in a column, as in some other models?

The sensitivity experiments nc00_w10 and nc00_w10k10 in group III in Table 1 were designed to provide qualitative answers to these questions. In the experiment nc00_w10, $w_{\text{act}}^{\ast }$, defined in Eq. (2), was artificially enhanced by a factor of 10 to affect aerosol activation but not turbulent mixing; in the experiment nc00_w10k10, both $w_{\text{act}}^{\ast }$ and the turbulent diffusivity coefficient, denoted 𝒦_mix, were enhanced by a factor of 10. To avoid causing substantial drifts of the mean climate, these changes were applied

• only in the DROPMIXNUC parameterization (i.e., with no code changes in the turbulent transport of heat and moisture parameterized with CLUBB),

• only in grid layers with pressure higher than 600 hPa, and

• only in grid columns where three or more layers in that altitude range had σ_w⩽ 0.1 m s⁻¹ and f_c⩾ 0.9.

The results are shown in Figs. 13 and 3.

Before discussing the results from these sensitivity experiments, it is worth pointing out an additional feature of the atmospheric conditions associated with ultra-low CDNCs. In Fig. 11g, the skewness of subgrid vertical velocity calculated by CLUBB is shown. Compared with the cases with higher CDNCs, ultra-low CDNCs are associated with significantly smaller skewness, especially in grid layers near and below ≈900 hPa. The smaller skewness of subgrid vertical velocity (Fig. 11g), weaker turbulence (Fig. 11f), larger cloud fractions (Fig. 10e), and lack of local cloud base in the grid boxes (Fig. 11d), together with the ParaView screenshots of Fig. 8 and the animation in the Video Supplement, give the impression that the ultra-low CDNCs in Arctic July often occur in lower-level stratus clouds with large horizontal and vertical extents and weak turbulence. The higher CDNC values are associated with larger skewnesses of subgrid vertical velocity (Fig. 11g), stronger turbulence (Fig. 11f), significant areas of local cloud base (Fig. 11d), and smaller cloud fractions (Fig. 10e), suggesting that these cases are more convective. These qualitative features are also seen in the other three regions, as shown in Figs. S10–S12. To further analyze the atmospheric environment, it would be useful to examine more physical quantities in the grid columns, for example, the surface fluxes, boundary layer stability, atmospheric conditions, and aerosol concentrations at the main cloud base. As mentioned earlier, the conditional sampling performed in this study was pointwise, hence we will defer the additional analyses to follow-on investigations.

6.2 Response of cloud droplet number

To present the responses of the simulated CDNCs to enhanced turbulence, we use here a metric ${\widehat{N}}_{\mathrm{l},\text{top}}$ from a protocol of the AeroCom community, namely, the in-cloud CDNC at the top of liquid water clouds. The top of liquid water clouds in a grid column is defined as the highest-in-altitude model layer with a grid box mean cloud liquid mass mixing ratio higher than 10⁻⁵ g kg⁻¹. As instructed by bullet one under “Q/A” at https://wiki.met.no/aerocom/indirect (last access: 31 August 2025), the time-averaged grid box mean CDNC and cloud fraction were diagnosed online and archived as model output. The annual mean ${\widehat{N}}_{\mathrm{l},\text{top}}$ was calculated by dividing the annual mean grid box mean liquid-cloud-top CDNC by the annual mean liquid-cloud-top cloud fraction.

Figure 13 (a–c) Annual mean in-cloud CDNC at the top of liquid water clouds, ${\widehat{N}}_{\mathrm{l},\text{top}}$, simulated under the PD emission scenario: (a) results from the nc00 configuration, (b) changes caused by enhancing the characteristic updraft velocity used in the aerosol activation calculation in the lower-tropospheric grid boxes with large cloud fractions and weak turbulence (nc00_w10 minus nc00), and (c) changes caused by enhancing both the characteristic updraft velocity and the diffusion coefficient for cloud droplet number (nc00_w10k10 minus nc00). (d–f) As (a–c), but showing the PD–PI differences in ${\widehat{N}}_{\mathrm{l},\text{top}}$. Masked out in white are locations with liquid cloud cover lower than 5 %. Details of the experimental design can be found in Sect. 6.1. Discussions of the results can be found in Sects. 6.2 and 6.3.

Figure 13a–c present the annual mean ${\widehat{N}}_{\mathrm{l},\text{top}}$ simulated under the PD emission scenario. Figure 13a shows the geographical distribution simulated by the control model, i.e., the nc00 configuration; Fig. 13b and c show the changes caused by enhancing $w_{\text{act}}^{\ast }$ alone or enhancing both $w_{\text{act}}^{\ast }$ and 𝒦_mix. Overall, the systematic changes in ${\widehat{N}}_{\mathrm{l},\text{top}}$ are positive, and Fig. 13c shows larger increases than Fig. 13b (note that the contour levels are logarithmic), suggesting that stronger aerosol activation and stronger droplet mixing can both help increase the simulated cloud droplet number. The largest systematic changes in ${\widehat{N}}_{\mathrm{l},\text{top}}$ are seen in the Arctic, SE China, the SE United States, and over the North Pacific storm track, matching some of the geographical regions identified in Sect. 4.2 with frequent occurrences of ultra-low CDNCs in EAMv2. This confirms that weak turbulence plays an important role in causing ultra-low CDNCs in the model. Furthermore, our analysis so far has led to the impression that the atmospheric conditions associated with ultra-low CDNCs exhibit features of low-level stratus clouds, while low-level stratus clouds are known to be commonplace in the aforementioned regions, suggesting that we have probably identified a cloud regime that is important for the occurrences of ultra-low CDNCs in the model.

Among the different geographical locations, the largest absolute changes in ${\widehat{N}}_{\mathrm{l},\text{top}}$ are found to occur in SE China and the SE United States. This can be understood from the equations and results presented earlier. In Sect. 2.3, it is explained that the characteristic updraft velocity, $w_{\text{act}}^{\ast }$, is used to calculate the fractions of interstitial aerosol particles being activated. From Eqs. (7) and (9), we see that the droplet number tendencies caused by the activation of interstitial APs in mode i, denoted here by 𝒰_i, can be expressed in an abstract form as

$$\begin{array}{}\text{(14)}& {\mathcal{U}}_i:={\left({\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{N}}_l}{\partial t}}\right)}_{\text{nuc},i}={\mathcal{G}}_i\left(w_{\text{act}}^{\ast }\right){\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}.\end{array}$$

The function 𝒢_i depends explicitly on $w_{\text{act}}^{\ast }$ (because of the activated fraction f_a,i), while ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}$ does not (at least to the first order). Hence, a change in this nucleation rate (denoted here by δ 𝒰_i) due to a change in $w_{\text{act}}^{\ast }$ (denoted here by $\mathit{\delta }{w}_{\text{act}}^{\ast }$) can be estimated using a Taylor expansion, namely,

$$\begin{array}{}\text{(15)}& \mathit{\delta }{\mathcal{U}}_i\approx \left({\displaystyle \frac{\partial {\mathcal{U}}_i}{\partial w_{\text{act}}^{\ast }}}\right)\mathit{\delta }{w}_{\text{act}}^{\ast }\approx \left({\displaystyle \frac{\partial {\mathcal{G}}_i}{\partial w_{\text{act}}^{\ast }}}\right)\left(\mathit{\delta }{w}_{\text{act}}^{\ast }\right){\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i},\end{array}$$

which is proportional to the interstitial AP concentration ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}$. Figure 12 shows that the ultra-low CDNCs in the polluted regions occur under higher ambient AP concentrations than in remote regions. Hence, given the same increases in $w_{\text{act}}^{\ast }$ (and the consequent increases in the activation fractions), the regions with higher AP concentrations are expected to feature larger changes in the number of nucleated droplets.

6.3 Responses of ERF_aer

Figure 13d–f present the PD–PI differences in ${\widehat{N}}_{\mathrm{l},\text{top}}$, again with Fig. 13d showing results from the nc00 configuration and Fig. 13e and f showing changes caused by increasing $w_{\text{act}}^{\ast }$ or both $w_{\text{act}}^{\ast }$ and 𝒦_mix. The key signals seen here are increases – not only in ${\widehat{N}}_{\mathrm{l},\text{top}}$ itself, as seen earlier in Fig. 13a–c, but also in the PD–PI differences of ${\widehat{N}}_{\mathrm{l},\text{top}}$, suggesting that there is a differential impact on ${\widehat{N}}_{\mathrm{l},\text{top}}$ under PD and PI emission scenarios. The reason for the differential impact can again be explained with Eq. (15): the PD emissions lead to higher ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{a},i}$ than the PI emissions and hence larger δ𝒰_i (nucleation rate increase) in response to the same $\mathit{\delta }{w}_{\text{act}}^{\ast }$. Keeping in mind the larger PI-to-PD increases in ${\widehat{N}}_{\mathrm{l},\text{top}}$ in the sensitivity experiments nc00_w10 and nc00_w10k10, it is not surprising to see in Fig. 3 that the global mean annual mean ERF_aer is also stronger in magnitude (more negative) in these experiments.

Recall that this study was motivated by the commonly observed model behavior that using a lower bound to remove extremely low CDNCs can help weaken the simulated global mean ERF_aer in aerosol-climate models. In the sensitivity experiments discussed in this section, CDNCs are increased in regions where ultra-low values are found in the nc00 configuration; hence, we expect ultra-low CDNCs to occur less frequently. The corresponding ERF_aer, however, strengthens rather than weakens in magnitude. We conducted additional experiments in which the assumed minimum value of the characteristic subgrid updraft velocity, i.e., w_act,min in Eq. (2), was changed. Those experiments exhibited the same qualitative behavior, namely, higher overall CDNCs are correlated with stronger (more negative) ERF_aer. This qualitative behavior has been observed in other models (e.g., Gantt et al., 2014), most recently in Ghosh et al. (2024). These results are undesirable from the ERF_aer perspective, but they do not contradict the weaker (less negative) ERF_aer seen in global models when lower bounds are applied to CDNC. At the locations and time steps where ultra-low CDNCs are simulated under both PD and PI emission scenarios, an imposed lower bound will artificially eliminate the PD–PI differences in CDNC and consequently eliminate their contribution to ERF_aer. We speculate that this is how the use of CDNC_min reduces ERF_aer in EAMv2 – and possibly in some other global models, as Mignot et al. (2021, Sect. 4.5) have discussed for the IPSL-CM6A-LR model.

It is also worth noting that our results, presented here, do not contradict the widely established notion of the cloud albedo susceptibility being weaker (less negative) in regimes with higher aerosol load and higher droplet number concentration. Figure 14 explains this point using a schematic in the spirit of Fig. 3 of Carslaw et al. (2013). The independent variable x and dependent variable y in Fig. 14 here can be interpreted as any of the pairs of variables depicted in Fig. 3b–d of Carslaw et al. (2013), although, for our discussion here, the most pertinent pairs are (I) CDNC as y with CCN concentration as x and (II) cloud albedo as y with anthropogenic aerosol emission as x. All three curves in Fig. 14 here, y₁, y₂, and y₃, have the same qualitative feature as in Fig. 3 of Carslaw et al. (2013), namely, their slopes ($\mathrm{d}y/\mathrm{d}x$) become less steep as x increases. The y₁ curve schematically represents the nc00 configuration of EAMv2 and the y₂ curve depicts our sensitivity experiment with enhanced $w_{\text{act}}^{\ast }$ (or $w_{\text{act}}^{\ast }$ as well as 𝒦_mix). Taking the variable pair (I) stated earlier as an example, Δx in our schematic is an emission-induced CCN concentration change; Δy₁ and Δy₂ are the resulting droplet number changes in the nc00 configuration and in the sensitivity experiments, respectively. Figure 13 shows that, with a larger $w_{\text{act}}^{\ast }$, both ${\widehat{N}}_{\mathrm{l},\text{top}}$ and its PD–PI difference increase systematically in many geographical regions. For Fig. 14, this means that the y₂ curve lies above the y₁ curve and has steeper slopes than the y₁ curve in the depicted range Δx. The y₃ curve symbolizes a hypothetical model configuration that simulates higher droplet concentrations than both nc00 and our sensitivity experiments with enhanced turbulence – so much higher that, for the CCN concentration range indicated by Δx, the simulation is much closer to the regime marked “droplet formation limited by supersaturation” in Fig. 3b of Carslaw et al. (2013), hence the y₃ curve has significantly gentler slopes that are similar to those found on the y₁ and y₂ curves at much larger values of x. The schematic illustrates that, if we focus on a single curve (which corresponds to a single model configuration with a specific choice of model formulation and parameters for cloud droplet nucleation), y is indeed less sensitive to x in regimes with higher x (and higher y, since the function is monotonic). However, different model configurations may appear as different curves in the schematic and larger y values do not necessarily correspond to smaller $\mathrm{d}y/\mathrm{d}x$ slopes across different curves. Hence, a model capable of simulating higher CDNCs does not necessarily produce weaker ERF_aer. This point can, in fact, be inferred from the uncertainty envelope diagram in Fig. 3e of Carslaw et al. (2013).

Figure 14 Schematic sketched in the spirit of Fig. 3 of Carslaw et al. (2013). The three curves, y₁, y₂, and y₃, symbolize three configurations of EAMv2 with differences in the parameters or model formulation used for cloud droplet nucleation in stratiform and shallow convective clouds. More details and discussions can be found in Sect. 6.3.

Results from both the literature and our experiments indicate that, while ultra-low CDNCs and strong ERF_aer are related, the relationship is complex. It is possible that the frequent occurrences of ultra-low CDNCs and relatively strong ERF_aer in EAMv2 have the same root causes but, very probably, there are also various more direct factors that affect CDNC and ERF_aer separately. Before the root causes are addressed and the more direct impactors are well understood, alleviation of one symptom might not automatically lead to mitigation of the other.

7 Summary and conclusions

This study aims at identifying where ultra-low CDNCs (i.e., in-cloud cloud droplet number concentrations lower than 10 cm⁻³) occur in the stratiform and shallow convective clouds simulated by E3SMv2 and which of the identified situations have the strongest impact on the simulated effective radiative effect of anthropogenic aerosols, ERF_aer. Process-level analyses are conducted to reveal characteristics of the cloud droplet formation, transport, and removal processes associated with impactful ultra-low CDNCs in stratiform and shallow convective clouds.

Simulations performed with present-day emissions show that ultra-low CDNCs occur most frequently over the mid- and high-latitude oceans in both hemispheres, with a general trend of higher frequencies at higher latitudes (Fig. 5a–c). Over the North Pacific storm track and in the Arctic region, among all the model time steps in a year and among all grid boxes in the middle and lower troposphere, about 20 %–40 % of the cases with cloud liquid presence are associated with ultra-low CDNCs (Fig. 6a and b). In polluted continental regions like SE China and the SE United States, this percentage is about 15 % or higher (Fig. 6c and d), despite the generally high aerosol concentrations (Fig. 12c and d).

Sensitivity experiments reveal that ultra-low CDNCs in the lower troposphere (with air pressure higher than 600 hPa) play the primary role in determining the sensitivity of E3SMv2's global mean ERF_aer to the lower bound of CDNC (CDNC_min = 10 cm⁻³, Fig. 3). The NH middle latitudes and the Arctic region show the largest absolute changes of annual mean ERF_aer when the lower bound is imposed (Fig. 4b), while the Arctic region features the largest relative changes in the annual mean ERF_aer (Fig. 4c). In the Arctic, the strongest seasonal responses of ERF_aer are found in boreal summer (Fig. 4f).

Ultra-low CDNCs are typically found in grid boxes with either very large (>0.9) or very small (<0.1) cloud fractions (Fig. 7); the former category is the larger contributor to the global mean ERF_aer's sensitivity to CDNC_min (Fig. 3). The ultra-low CDNCs associated with large cloud fractions are typically found in contiguous grid boxes in the main body of large cloudy regions (Figs. 8 and 9 and the Video Supplement). These cases feature strong water vapor condensation (Fig. 10g, red line) and weak turbulence (Fig. 11f, red line), while the statistical distributions of subgrid vertical velocity feature relatively small skewness (Fig. 11g, red line), suggesting that the clouds are likely to be low-level stratus.

A comparison of the cloud droplet number budgets associated with different CDNC ranges suggests that the ultra-low CDNCs are caused by weak sources of droplet number (Fig. 10c), especially lack of aerosol activation (Fig. 11a and b), rather than strong sinks of droplet number (Figs. 10b and d and 11c).

CCN number concentrations were diagnosed by calculating the number concentrations of aerosol particles of all sizes and compositions in the grid box that would be activated at 0.1 % supersaturation. In March over SE China, the composite mean CCN concentrations associated with ultra-low CDNCs are of the order of several hundred per cm³ (Fig. 12c), suggesting that the simulated cloud droplet nucleation process (rather than the aerosol particle number concentration) plays the primary role in causing the ultra-low CDNCs in that region and time of year. In the selected focus months over the North Pacific storm track, the Arctic region, and the SE United States, the composite mean CCN number concentrations associated with ultra-low CDNCs range from close to 20 cm⁻³ to about 120 cm⁻³ in the lower troposphere and the values are significantly lower than the composite mean CCN concentrations associated with higher CDNCs (Fig. 12a, b, and d), suggesting that the simulated cloud droplet nucleation process and aerosol concentrations are both possible contributors to the ultra-low CDNCs.

Sensitivity experiments show that, in lower-tropospheric grid boxes with large cloud fractions and weak turbulence, boosting the characteristic updraft velocity used in the aerosol activation parameterization or enhancing the turbulent mixing of cloud droplet number can increase the simulated CDNCs (Fig. 13b and c) in regions where ultra-low CDNCs are frequently found in EAMv2 (Fig. 5a–c) and where the use of CDNC_min significantly affects ERF_aer (Fig. 4b), suggesting that the identified atmospheric conditions are relevant. However, those increased CDNCs are accompanied by increased (rather than decreased) PD–PI differences in CDNC (Fig. 13e and f), as well as larger (rather than smaller) magnitudes of the global mean ERF_aer (Fig. 3). This result is qualitatively consistent with findings from earlier studies, such as Gantt et al. (2014), and most recently from Ghosh et al. (2024). Our interpretation of this undesirable correlation between changes in CDNC and changes in ERF_aer is that, while the frequent occurrences of ultra-low CDNCs and relatively strong ERF_aer in EAMv2 may have the same root causes, there are probably more direct factors that affect CDNC and ERF_aer separately. Here, we refrain from speculating what the root causes and more direct factors might be, noting that the E3SM model developers are actively working in this area (Shan et al., 2024a, b; Burrows et al., 2024). Nevertheless, our results, presented in this paper, clearly suggest that mid- and high-latitude low-level stratus occurring under weak turbulence is a cloud regime worth further investigating.