Magnitude and timescale of liquid water path adjustments to cloud droplet number concentration perturbations for nocturnal non-precipitating marine stratocumulus

Cloud liquid water path (L) adjusts to perturbations in cloud droplet number concentration (N) over time. We explore the magnitude and timescale of this adjustment in nocturnal non-precipitating marine stratocumuli using large eddy simulations of baseline conditions and aerosol seeding experiments for 22 meteorological conditions. The results confirm that the L adjustment (δL) slope (k) is more negative for simulation pairs with relatively low N and less negative for high N. Overall, k is unlikely to be lower than −0.4 within 24 h of seeding start, meaning the L adjustment is unlikely to fully offset the brightening due to the Twomey effect. After seeding, the δL becomes increasingly negative, which can be characterized by an exponential convergence. This evolution is governed by a short timescale around 5 h and lasts for around 8–12 h. It is driven by the feedback between entrainment, L, and boundary layer (BL) turbulence. Other processes, including radiation, surface fluxes, and subsidence, respond to the seeding weakly. This short timescale is insensitive to the amount of seeding, making the evolution of δL and some other deviations similar for different seeding amounts after appropriate scaling. The timescale of k evolution is closely related to the δL timescale and hence also short, while it could also be affected by the δN evolution. The results are most relevant to conditions where seeding is applied to a large area of marine stratocumulus in well-mixed and overcast BL where shear is not a primary source of turbulence.

1 Introduction

Marine cloud brightening (MCB) has been proposed as a climate intervention strategy to mitigate global warming by taking advantage of the Twomey effect (Twomey, 1974, 1977): injecting aerosol particles into marine stratocumuli (henceforth “seeding”) increases cloud droplet number concentration (denoted by N) and reduces cloud droplet size; with cloud water amount unchanged, this perturbation increases the cloud albedo to reflect more solar radiation back to space to cool the Earth (Latham, 1990; Feingold et al., 2024). The initial brightening occurs within a short space of time around 10–15 min as the seeded particles are transported vertically through the boundary layer (BL). Afterwards, the macroscopic properties of these clouds, including liquid water path (denoted with L) and cloud fraction, adjust to the changes in N (Albrecht, 1989; Ackerman et al., 2004; Bretherton et al., 2007), further modifying the albedo of a cloudy scene. The sign, the magnitude, and the timescale of these adjustments may affect the effectiveness of MCB.

In this study, we focus on the adjustment of L to a perturbation in N (abbreviated as “L adjustment”) that is initially caused by seeding. This adjustment can be characterized by the ratio

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }k={\displaystyle \frac{\mathit{\delta }l}{\mathit{\delta }n}},\end{array}$$

where l=ln L, n=ln N, and δ indicates the difference between seeded and unseeded conditions. We refer to k as the adjustment slope because it is the slope of the line segment connecting cloud states (N,L) with and without seeding in the N–L plane on a log scale. It can be used to characterize the sensitivity of l to n required to assess the susceptibility of cloud albedo to N under idealized conditions (Platnick and Twomey, 1994; Bellouin et al., 2020). A positive k indicates that the L adjustment further brightens a cloudy scene in addition to the Twomey effect, while a negative k indicates the opposite. In particular, a k of −0.4 indicates that the L adjustment exactly offsets the Twomey effect (note that this value of k sets the right-hand side of Eq. 13 in Bellouin et al., 2020, to zero, indicating no change in cloud albedo).

For non-precipitating marine stratocumuli, k is likely negative but its magnitude is still uncertain. The negative sign comes from a few mechanisms. An initial increase in N leads to smaller droplets. They evaporate faster during the mixing between the cloud and the free-troposphere (FT) air (Wang et al., 2003); they sediment slower, allowing more liquid near the cloud top to evaporate (Bretherton et al., 2007). The reduction of the drop size also suppresses the cloud-base precipitation, which may strengthen the turbulence (Lu and Seinfeld, 2005; Caldwell et al., 2005; Wood, 2007; Sandu et al., 2008). Smaller droplets may also enhance cloud-top longwave radiative cooling (Garrett et al., 2002; Petters et al., 2012; Igel, 2024). All these mechanisms suggest that an increase in N enhances the cloud-top entrainment, which leads to more warming and drying of the BL and reduces L, although there is still debate on the dominant mechanism behind this enhancement (Igel, 2024). During the daytime, higher N causes more absorption of solar radiation (Stephens, 1978; Boers and Mitchell, 1994; Petters et al., 2012), also contributing to a negative k, although the weaker absorption by lowered L may limit this effect (Sandu et al., 2008; Zhang et al., 2024). The uncertainty in our understanding of the magnitude of k can be seen in the wide range of values reported in previous studies, which used different methods and focused on different conditions (see, e.g., a summary compiled in Fig. 1 and Table S1 in the Supplementary Material of Glassmeier et al., 2021). The correlation between environmental conditions and k is also under debate (e.g., Chun et al., 2023).

Several recent works (e.g., Gryspeerdt et al., 2021; Glassmeier et al., 2021) brought attention to the timescale for L to adjust. Since

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\dot{k}={\displaystyle \frac{\mathit{\delta }\dot{l}}{\mathit{\delta }n}}-k{\displaystyle \frac{\mathit{\delta }\dot{n}}{\mathit{\delta }n}},\end{array}$$

where the dot above a variable indicates the time derivative ($\mathrm{d}/\mathrm{d}t$), the timescales for both L adjustment and N evolution contribute to the timescale for k (Gryspeerdt et al., 2022). Glassmeier et al. (2021) reported that k approaches a “steady-state” slope (k_∞) around −0.64 with an adjustment timescale around 20 h based on the analysis of a large eddy simulation (LES) ensemble of more than 100 nocturnal marine stratocumuli. Other studies have reported shorter adjustment timescales (e.g., Rahu et al., 2022; Prabhakaran et al., 2023). Reconciling different estimates of adjustment timescales may require connecting them to various other timescales previously discovered in marine stratocumulus-topped BLs (STBLs), e.g., the inversion adjustment timescale (Schubert et al., 1979), the thermodynamic timescale (Schubert et al., 1979; Bretherton et al., 2010), and a short timescale associated with the strong feedback between entrainment velocity, L, and BL turbulence (Zhu et al., 2005; Bretherton and Blossey, 2014) as shown by Jones et al. (2014).

One less investigated aspect of the L adjustment in the non-precipitating regime is its dependence on N. The effects of several aforementioned mechanisms for entrainment enhancement by increasing N should saturate at high N where a further increase in N does not significantly affect the drop size, suggesting that the negative k caused by these mechanisms at relatively low N should become less negative at high N. Results from some previous studies (Lu and Seinfeld, 2005; Chen et al., 2011) support this expectation. Further narrowing down the uncertainty in the N dependence of k is valuable because it may reveal optimal MCB seeding strategies and provide an assessment of the integrated brightening effects as the seeded clouds lose N through dilution (e.g., by mixing with FT air).

In this study, we use LES to investigate the magnitude and timescale of the L adjustment for nocturnal non-precipitating marine stratocumuli. Even though applications like MCB intrinsically involve cloud evolution during the day, the nighttime evolution is still important because the high-N conditions due to previous seeding may continue into the nighttime, and the nighttime evolutions of the clouds and the BL shape the evolution over the course of the following day (Sandu et al., 2008; Chun et al., 2023; Zhang et al., 2024). This choice of scope also reduces the processes involved and simplifies the problem.

2 Method

2.1 Model and shared simulation configurations

LESs are performed using the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall, 2003), version 6.10.10. The choices of specific schemes (e.g., numerical schemes for the dynamics, physical parameterizations) are identical to those described in the first two paragraphs of Sect. 2 in Chen et al. (2024).

Simulation configurations largely follow the configurations used in Glassmeier et al. (2021). The simulation domain is $\mathrm{48}\mathrm{km}\times \mathrm{48}\mathrm{km}\times \mathrm{2}\mathrm{km}$ in the x, y, and z dimensions with 200 m horizontal and 10 m vertical grid spacings. Periodic lateral boundary conditions and a damping layer from 1.5 km to domain top are used. The subsidence profile follows

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{w}_{\mathrm{s}}=\left\{\begin{array}{ll}-Dz,& z<\mathrm{1600}\mathrm{m}\\ \mathrm{0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},& z\ge \mathrm{1600}\mathrm{m},\end{array}\right.\end{array}$$

where the divergence $D=\mathrm{3.75}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻¹ is based on the DYCOMS-II RF02 case (Wyant et al., 2007; Ackerman et al., 2009). The large-scale wind speed is zero at all levels so that there is no mean shear to drive turbulence. A constant surface aerosol flux of 70 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ is prescribed to be consistent with Glassmeier et al. (2021). The time step is 1 s and the radiative scheme is called once every 10 s.

The main differences from Glassmeier et al. (2021) lie in other lower boundary conditions. We calculate surface fluxes interactively, instead of prescribing constant values based on DYCOMS-II RF02. A wind speed of 7 m s⁻¹, based on ERA5 climatology (Hersbach et al., 2020), is added to the surface local wind fluctuation when calculating surface sensible and latent heat fluxes. See details in Sect. 2 and Appendix A in Chen et al. (2024). The sea surface temperature (SST) is case-dependent and fixed at 0.5 K warmer than the initial surface air temperature, instead of being fixed at one value for all simulations.

2.2 Case definition and experiment design

We arbitrarily select 22 non-precipitating cases from Glassmeier et al. (2021). As shown in Fig. S1 in the Supplement, they span a wide range in the N–L plane and the plane of L and inversion base height (z_i). Here and throughout the paper, N is characterized by the mean cloud droplet number concentration for cloudy columns with hydrometeor optical depth greater than 1 and L by the domain mean LWP, both following Glassmeier et al. (2021). In Glassmeier et al. (2021), a case is defined by its initial thermodynamic and aerosol profiles that are controlled by six parameters: initial boundary layer (BL) depth (h_mix) randomly drawn from 500 to 1300 m, BL liquid water potential temperature (θ_l) from 284 to 294 K, BL total water mixing ratio (q_t) from 6.5 to 10.5 g kg⁻¹, θ_l jump across the inversion base (Δθ_l) from 6 to 10 K, q_t jump (Δq_t) from −10 to −6 g kg⁻¹, and aerosol mixing ratio (N_a) from 30 to 500 mg⁻¹. As shown in Fig. S2, our 22 cases also cover the pair-wise spaces of the first five parameters, which are hereafter collectively referred to as the meteorological condition (MC). Note that the ensemble in Glassmeier et al. (2021) is characterized by a relatively dry free troposphere (FT) with FT q_t between 0.2 and 2.8 g kg⁻¹ (see the BL q_t–Δq_t panel in Fig. S2). Our MCs are the same in this respect.

We use the MCs from these 22 cases to set up the initial θ_l and q_t profiles but configure the aerosol as follows. For each MC, we first find a BASE run by setting the initial N_a throughout the domain to a low value and perform a 36 h simulation so that the trajectory of the simulation in the N–L plane is close to the approximate threshold for precipitation (defined as the line that corresponds to a characteristic cloud-top mean drop radius of 12 µm) but the simulated cloud remains non-precipitating (defined as a cloud-base precipitation rate of less than 0.5 mm d⁻¹; Wood, 2012). Then, we seed the BASE run at 12 h after the beginning of the simulation by uniformly increasing the total number mixing ratio (i.e., the sum of aerosol, cloud droplet, and raindrop number mixing ratios) from the surface to 400 m for 30 min at constant rates estimated to achieve seven N targets: 95, 125, 165, 225, 300, 400, 550 cm⁻³. A seeded run for a given MC is only performed if the N target is greater than N in the BASE run at the end of its 36 h simulation. The seeded runs are then simulated for 24 h. Due to various processes affecting N, the actual N values for seeded runs aimed at these seven N targets reach about 120, 135, 165, 200, 250, 310, and 400 cm⁻³ at 36 h, still spanning a relatively wide range. We refer to seeded runs by their final N as N120, N135, N165, N200, N250, N310, and N400, respectively.

The simulation set described above is referred to as MAIN. Additional sensitivity runs are performed and will be introduced as needed.

2.3 Budget analysis

2.3.1 L budget

The L budget is diagnosed based on mixed-layer theory (MLT; Lilly, 1968; Wood, 2007; van der Dussen et al., 2014; Ghonima et al., 2015; Hoffmann et al., 2020) because the BLs in all simulations are close to well-mixed and overcast (Fig. S3). We start with

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }L={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{\Gamma }}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle (z_{\mathrm{i}}-z_{\text{cb}}{)}^{\mathrm{2}},\end{array}$$

where Γ_l is the rate of change of the liquid water mixing ratio (q_l) with height in an adiabatic cloud, 〈ρ₀〉 is the cloud layer mean air density, and z_cb is the domain mean cloud base. Then,

$$\begin{array}{ll}{\displaystyle }\dot{L}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\dot{\mathrm{\Gamma }}}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle \right)(z_{\mathrm{i}}-z_{\text{cb}}{)}^{\mathrm{2}}+{\mathrm{\Gamma }}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle (z_{\mathrm{i}}-z_{\text{cb}})\left[{\dot{z}}_{\mathrm{i}}-{\dot{z}}_{\text{cb}}\right]\\ {\displaystyle }& {\displaystyle }\approx {\mathrm{\Gamma }}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle (z_{\mathrm{i}}-z_{\text{cb}})\left[{\dot{z}}_{\mathrm{i}}-{\dot{z}}_{\text{cb}}\right]\\ \text{(5)}& {\displaystyle }& {\displaystyle }\approx {\mathrm{\Gamma }}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle (z_{\mathrm{i}}-z_{\text{cb}})\left[{\dot{z}}_{\mathrm{i}}-\left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}\dot{\langle q_{\mathrm{t}}\rangle }\right)\right],\end{array}$$

where 〈θ_l〉 and 〈q_t〉 are the BL mean θ_l and q_t, and $\partial z_{\text{cb}}/\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle$ and $\partial z_{\text{cb}}/\partial \langle q_{\mathrm{t}}\rangle$ are the sensitivities of z_cb to BL warming/cooling and drying/moistening. Although Γ_l, 〈ρ₀〉, $\partial z_{\text{cb}}/\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle$, and $\partial z_{\text{cb}}/\partial \langle q_{\mathrm{t}}\rangle$ also change slightly with time, we use their time-dependent values in Eq. (5) but ignore terms containing their temporal rates of change when expanding $\dot{L}$. The contributions of various physical processes to $\dot{L}$ are then diagnosed via their effects on $\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }$, $\dot{\langle q_{\mathrm{t}}\rangle }$, and ${\dot{z}}_{\mathrm{i}}$:

$$\begin{array}{ll}{\displaystyle }{\dot{L}}_P=& {\displaystyle }{\mathrm{\Gamma }}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle (z_{\mathrm{i}}-z_{\text{cb}})\\ \text{(6)}& {\displaystyle }& {\displaystyle }\times \left[{\dot{z}}_{\mathrm{i},P}-\left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}{\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}_P+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}{\dot{\langle q_{\mathrm{t}}\rangle }}_P\right)\right],\end{array}$$

where tendency terms with a subscript P indicate the contributions of P, which refers to cloud-top entrainment (ENTR), radiation (RAD), subsidence (SUBS), surface fluxes (SURF), or precipitation (PRCP).

Regarding $\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }$ and $\dot{\langle q_{\mathrm{t}}\rangle }$, we express the contributions of a process in the form of the BL flux divergences of θ_l and q_t ($\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,P}$ and $\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,P}$) evenly distributed over the entire BL depth, assuming that the BL is well-mixed,

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}_P={\displaystyle \frac{\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,P}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}}},\dot{\langle q_{\mathrm{t}}\rangle }={\displaystyle \frac{\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,P}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}}}.\end{array}$$

For RAD, $\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,\text{RAD}}$ is the radiative heating rate integrated from the surface to z_i. For SURF, the corresponding fluxes (i.e., surface sensible and latent heat fluxes) are only defined at the surface, but we still express its contributions in terms of flux divergences with the corresponding fluxes at z_i defined to be zero. Similarly, the flux divergences of PRCP are based on the surface precipitation rate and a zero precipitation flux at z_i. SUBS is assumed to have no contributions to either $\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }$ or $\dot{\langle q_{\mathrm{t}}\rangle }$ because the subsidence as prescribed in Eq. (3) only stretches the vertical profiles and does not move any air parcel across the z_i. For ENTR, we simply attribute the differences between actual $\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }$ and $\dot{\langle q_{\mathrm{t}}\rangle }$ and the sums of contributions of RAD, SURF, and PRCP to ENTR, closing the 〈θ_l〉 and 〈q_t〉 budgets without residuals (see Sect. 4.1 in Chen et al., 2024, for details). We still diagnose $\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,\text{ENTR}}$ and $\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,\text{ENTR}}$ from ${\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}_{\text{ENTR}}$ and ${\dot{\langle q_{\mathrm{t}}\rangle }}_{\text{ENTR}}$ following Eq. (7), which will be used later in the paper.

Regarding ${\dot{z}}_{\mathrm{i}}$, the SUBS term is simply the w_s(z_i) calculated using Eq. (3) and the ENTR term is the entrainment velocity:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{w}_{\mathrm{e}}={\dot{z}}_{\mathrm{i}}-w_{\mathrm{s}}\left(z_{\mathrm{i}}\right).\end{array}$$

Finally, a residual term (RES) is calculated from the difference between the actual $\dot{L}$ and the sum of diagnosed terms to close the L budget, even though it is not necessary for $\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }$, $\dot{\langle q_{\mathrm{t}}\rangle }$, or ${\dot{z}}_{\mathrm{i}}$.

2.3.2 L budget in terms of BL flux divergence

We define

$$\begin{array}{ll}{\displaystyle }{\mathcal{F}}_P=& {\displaystyle }\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}{\dot{z}}_{\mathrm{i},P}\\ \text{(9)}& {\displaystyle }& {\displaystyle }-\left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,P}+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,P}\right).\end{array}$$

Obviously, the sum of ℱ_P over ENTR, RAD, SUBS, SURF, and PRCP is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathcal{F}=\sum _P{\mathcal{F}}_P\\ {\displaystyle }& {\displaystyle }=\sum _P\left(\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}{\dot{z}}_{\mathrm{i},P}-\left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,P}+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,P}\right)\right)\\ \text{(10)}& {\displaystyle }& {\displaystyle }=\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}\left[{\dot{z}}_{\mathrm{i}}-\left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}\dot{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle }+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}\dot{\langle q_{\mathrm{t}}\rangle }\right)\right].\end{array}$$

Compared with Eq. (5), ℱ is essentially ${\dot{z}}_{\mathrm{i}}-{\dot{z}}_{\text{cb}}$ multiplied by 〈ρ₀〉_BLz_i. Hereafter we refer to ℱ as “the flux divergence for $\dot{L}$” to emphasize that we formulate it and its components by processes in terms of flux divergences $\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,P}$ and $\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,P}$.

Then, $\dot{L}$ can be written as

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\dot{L}=p{\mathit{\zeta }}_{\mathrm{c}}\mathcal{F},\end{array}$$

where

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }p={\mathrm{\Gamma }}_{\mathrm{l}}\langle {\mathit{\rho }}_{\mathrm{0}}\rangle /\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}\end{array}$$

is a prefactor and

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\mathit{\zeta }}_{\mathrm{c}}=\mathrm{1}-z_{\text{cb}}/z_{\mathrm{i}}\end{array}$$

is the normalized cloud depth.

2.3.3 k evolution

We define

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{\dot{k}}_L={\displaystyle \frac{\mathit{\delta }\dot{l}}{\mathit{\delta }n}}=k{\displaystyle \frac{\mathit{\delta }\dot{l}}{\mathit{\delta }l}}\end{array}$$

and

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{\dot{k}}_N=-k{\displaystyle \frac{\mathit{\delta }\dot{n}}{\mathit{\delta }n}}.\end{array}$$

Then, $\dot{k}$ is the sum of ${\dot{k}}_L$ and ${\dot{k}}_N$ (Eq. 2). And ${\dot{k}}_L$ can be further connected to the L budget terms:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{\dot{k}}_L=\sum _P{\dot{k}}_P=\sum _P{\displaystyle \frac{\mathit{\delta }{\dot{l}}_P}{\mathit{\delta }n}}=\sum _P{\displaystyle \frac{\mathit{\delta }\left({\dot{L}}_P/L\right)}{\mathit{\delta }n}},\end{array}$$

where P refers to ENTR, RAD, SUBS, SURF, PRCP, or RES. In this paper, we do not further decompose ${\dot{k}}_N$.

3 Results

In this section, we show results based on the simulations described in Sect. 2.2 (the MAIN set). We first present an overview of our simulations, then show the L and k budgets, and end with results about timescales of L, δL, and k. In this section, all times refer to the time from the moment when seeding starts, t_seeding, which is 12 h from the beginning of the simulations.

3.1 Overview

We present an overview starting with the cloud states (N, L) and the L adjustment slopes k in two 1 h windows (Fig. 1). Between 11 and 12 h, the 1 h median L tends to be lower for higher N for all MCs (Fig. 1), consistent with the expected negative L adjustment slope for non-precipitating stratocumuli.

Figure 1 Overview of 1 h medians of cloud states (N, L) and L adjustment slopes k between 11 and 12 h (a, b) and between 23 and 24 h (c, d) in MAIN. In panels (a) and (c), each curve connects simulations for one MC and there are at least six symbols along each curve (BASE, N165, N200, N250, N310, and N400). If present, symbols between BASE and N165 indicate N120 and N135. The combinations of the color and the shape of the symbol uniquely identify the 22 MCs, matching those in Figs. S1 and S2. The dash-dotted black lines are reference lines showing k at 12 and 24 h assuming it exponentially converges to −0.64 with a timescale of 20 h, based on Glassmeier et al. (2021). The dashed black line is the precipitation line defined as the line that corresponds to a characteristic cloud-top mean drop radius of 12 µm. Panels (b) and (d) show distributions of 1 h median k and 1 h median mid-point N for various groups of simulation pairs. Blue symbols and lines are based on k between BASE and seeded runs; red symbols are based on k for line segments in panels (a) and (c). See annotations in matching colors.

The distributions of k during the same time period are shown in Fig. 1b. Since BASE, N165, N200, N250, N310, and N400 are always simulated for all 22 MCs, we divide the slopes between BASE and various seeded runs into six groups: BASE–N165, BASE–N200, BASE–N250, BASE–N310, BASE–N400, and BASE–other seeded runs. Each of the first five groups includes, obviously, 22 pairs of simulations, while the last one also includes 22 pairs (i.e., 4 BASE–N120 pairs and 18 BASE–N135 pairs). As shown with blue symbols in Fig. 1b, the 1 h median k between the simulation pairs becomes less negative for larger N. The medians for the six groups increase from −0.38 for BASE–other seeded runs, which is more negative than predicted by Glassmeier et al. (2021), to −0.19 for BASE–N400.

This N dependence of k can also be seen in k for all pairs of simulations with adjacent N or, in other words, all line segments in Fig. 1a. We divide these line segments into five groups: N165–N200, N200–N250, N250–N310, and N310–N400, each including 22 line segments, and all remaining 44 line segments (i.e., 4 BASE–N120 segments, 14 BASE–N135 segments, 4 BASE–N165 segments, 4 N120–N135 segments, and 18 N135–N165 segments). The distributions of 1 h median k for line segments in five groups are shown with red symbols in Fig. 1b. The slopes are also less negative for pairs with larger N with medians ranging from −0.29 to −0.13. Even though the k values for high-N segments are between two seeded runs, not a BASE run with relatively high N and a seeded run from this BASE run, we can generally claim that the k values for high-N pairs are less negative than for low-N pairs. See Appendix A for more details on this point.

Between 23 and 24 h, the overall negative L adjustment slope does not change much from between 11 and 12 h (cf. Fig. 1c and a). The 1 h median k values between BASE and seeded runs (blue symbols) and for all line segments (red symbols) become more negative by about 0.06 (cf. Fig. 1d and b), but k values for most pairs are less negative than predicted by Glassmeier et al. (2021). The distributions of k in all groups are broader than 12 h earlier, partially because L for simulations within each MC fluctuates more (see, e.g., polylines with red × symbols and purple triangles in Fig. 1c).

In all these results, k is unlikely to reach −0.4 to fully compensate for the Twomey effect.

Next, we present the time series of L, N, and k, focusing on the 22-MC composites (i.e., averages) for four aerosol configurations: BASE, N165, N250, and N400 (Fig. 2). For BASE, L increases from 71 to 92 g m⁻² from 0 to 24 h since t_seeding. N increases from 107 to 125 cm⁻³, meaning the surface aerosol flux, which is the only aerosol source in BASE, dominates the net trend in $\dot{N}$. For seeded runs, L and N deviate from those for BASE. L for N165 still increases, but at a slower rate; L for N250 and N400 decreases first and then increases at even slower rates. N peaks soon after 1 h, overshooting the N target (165, 300, and 550 cm⁻³ for N165, N250, and N400) by a few percent, then decreases relatively quickly for about 2 h. This behavior is driven by the activation of the seeded aerosol and some initial adjustment of the boundary layer (BL) to the seeding. After about 5 h, N decreases steadily towards the end of the simulations. For N400, this indicates the dominance of the dilution by the mixing between BL and FT air, given the negligible collision–coalescence at high N; for N165, multiple processes balance each other and N does not change much. The k values between BASE and seeded runs and between seeded runs (blue and red curves in Fig. 2c) become negative after the seeding starts, initially quickly and then more slowly. At 24 h, L, N, and k for most aerosol configurations are not steady.

Figure 2 Time series of (a) L, (b) N, and (c) k averaged across simulations for BASE, N165, N250, and N400.

3.2 L budget

In this subsection, we examine the time series of the L budget, again focusing on the 22-MC composites for BASE, N165, N250, and N400. Since the budget terms are noisier than L, the results are smoothed with a 4 h running average and plotted at the end of the 4 h window. Recall that time is relative to t_seeding; i.e., curves before 4 h are affected by results before the seeding, which are the same for BASE, N165, N250, and N400. Due to the smoothing, the data shown for the first few hours since t_seeding are not suitable for directly estimating timescales, which will be investigated with unsmoothed data in Sect. 3.4, but the main features described in this subsection are valid.

We start with the time series of the flux divergence for $\dot{L}$ (ℱ) in Fig. 3a. The variations in time and between BASE and seeded runs are much smaller than the range of values spanned by different terms. To see details, we also show all combinations of budget terms and aerosol configurations as the differences from the values at 0 h in Fig. 3b and c. As expected, the contributions of ENTR and SUBS are negative (thinning cloud layer) and those of RAD and SURF are positive (deepening cloud layer). Since all simulated clouds are non-precipitating, the PRCP term is calculated for a precise RES term but omitted from the plot. The RES is close to 0, indicating a good closure of the budget (Fig. 3b).

Figure 3 The 4 h smoothed time series of variables related to the flux divergence for $\dot{L}$ (ℱ) and the normalized cloud depth (${\mathit{\zeta }}_{\mathrm{c}}=\mathrm{1}-z_{\text{cb}}/z_{\mathrm{i}}$), averaged across simulations for BASE, N165, N250, and N400: (a) actual ℱ and its budget terms minus values at t_seeding; (b) actual ℱ and its ENTR and RES terms minus values at t_seeding; (c) RAD, SUBS, SURF, and the sum of these three non-ENTR terms for ℱ; and (d) ζ_c.

For seeded runs, ℱ for N400 becomes less positive than for BASE after the seeding starts, then turns negative, and then slowly returns to the ℱ for BASE (Fig. 3a). The response in ℱ to seeding is dominated by the ENTR term (Fig. 3b). This is consistent with our understanding that the entrainment is enhanced by an increase in N. The RAD term is slightly weaker (less positive) and the SURF and the SUBS terms are slightly stronger (more positive and more negative) in seeded runs, but the sum of all three non-ENTR terms is slightly weaker (less positive). Even though the 4 h running average smooths the time series, it is clear that the response in SURF (Fig. 3c) is delayed compared with the responses in ENTR and RAD. N165 and N250 show similar behavior with weaker deviation from BASE.

The normalized cloud depth ζ_c for seeded runs become less than BASE (Fig. 3d), mainly because the cloud layers become thinner than BASE and not because of the faster growth of z_i due to enhanced entrainment (shown later).

We take a closer look at time series relevant to the entrainment since it dominates the response in ℱ. The entrainment velocity w_e dominates the response in ${\dot{z}}_{\mathrm{i}}$ to seeding (Fig. 4a). Initially, the seeded runs entrain much more strongly than BASE, causing faster growth of z_i; after a few hours, this difference starts to narrow. This behavior is similar to the results from Prabhakaran et al. (2023) (see their Fig. 3g). Since the Δθ_l and Δq_t in BASE only change slightly over time (Fig. S4a) and the Δθ_l and Δq_t respond weakly to seeding (Fig. S4b), the behavior of w_e dominates the response in ℱ_ENTR to seeding. At the end of the simulation, w_e in N400 is faster than in BASE by less than 0.2 mm s⁻¹, while the subsidence velocity w_s is more negative by less than 0.1 mm s⁻¹. As a result, z_i only grows marginally faster in seeded runs at 24 h since t_seeding.

Figure 4 Same as Fig. 3 but for variables related to the entrainment velocity (w_e): (a) ${\dot{z}}_{\mathrm{i}}$, w_e, and w_s minus values at t_seeding and (b) components in w_e parameterization $w_{\mathrm{e}}=(A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}})w_{\ast }^{\mathrm{3}}/z_{\mathrm{i}}$ divided by values at t_seeding. See text around Eqs. (17) and (18) for definitions of symbols.

We further break down w_e following

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{w}_{\mathrm{e}}=\left({\displaystyle \frac{A}{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}}}\right){\displaystyle \frac{w_{\ast }^{\mathrm{3}}}{z_{\mathrm{i}}}},\end{array}$$

where A is an entrainment efficiency, Δθ_v is the virtual temperature jump across z_i, and w_∗ is the convective velocity scale defined as

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{w}_{\ast }={\left({\displaystyle \frac{\mathrm{2.5}g}{T_{\mathrm{0}}}}\underset{\mathrm{0}}{\overset{z_{\mathrm{i}}}{\int }}\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}\mathrm{d}z\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where g is gravitational acceleration and T₀ is a temperature scale. We use $w_{\ast }^{\mathrm{3}}$ and z_i diagnosed by SAM and calculate $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$ as one term based on Eq. (17). The evolution of $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$ and its response to seeding are dominated by A, given the weak response in the Δθ_v inferred from the behavior of Δθ_l and Δq_t (Fig. S4). As shown in Fig. 4b, $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$ in BASE increases over time as both L and N in BASE increase, consistent with previous works (Nicholls and Turton, 1986; Turton and Nicholls, 1987; Bretherton et al., 2007); $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$ in seeded runs increases with the amount of seeding, meaning the enhancement by increasing N dominates over decreasing L. The evolution of $w_{\ast }^{\mathrm{3}}$ in BASE and its response to seeding are consistent with the positive correlation between the buoyancy flux and L (Bretherton and Wyant, 1997). The responses in $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$ and $w_{\ast }^{\mathrm{3}}$ to seeding are consistent with Chun et al. (2023). Different from w_e, the magnitudes of the responses in both $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$ and $w_{\ast }^{\mathrm{3}}$ to seeding remain relatively large throughout the simulations.

Figure 5 Same as Fig. 3 but for the (a) L budget and (b) l budget.

Figure 5a shows the time series of L budget terms. Since the prefactor p in Eq. (11) evolves slowly in BASE and does not respond much to the seeding (not shown), the response in the L budget to seeding is driven by ℱ and ζ_c. The much more negative ENTR term for ℱ in the seeded runs translates to a slightly more negative ENTR term for $\dot{L}$ (cf. Figs. 3a and 5a). Later in the simulation, the contributions of ENTR, RAD, SUBS, and SURF are all weaker in the seeded runs due to smaller ζ_c. Figure 5b shows the time series of the l budget terms. They are simply diagnosed as $\dot{l}=\dot{L}/L$ but still worth showing because $\dot{l}$ directly affects the evolution of k (Eq. 14). With seeding, the contributions of all processes become stronger, opposite to the responses in the L budget terms, because of the different ζ_c dependence between $\dot{l}$ ($\propto {\mathit{\zeta }}_{\mathrm{c}}^{-\mathrm{1}}\mathcal{F}$) and $\dot{L}$ (∝ζ_cℱ).

To summarize, entrainment dominates the response in ℱ to seeding. The response in ℱ_ENTR is dominated by the response in w_e, which peaks quickly and then decays towards the end of the simulation as a result of the balance between the persisting responses in A and $w_{\ast }^{\mathrm{3}}$. The responses in ℱ and ζ_c shape the evolution of $\dot{L}$ and $\dot{l}$.

3.3 k budget

In this subsection, we examine the evolution of k for two pairs of aerosol configurations: one between BASE and N165 and the other one between N250 and N400. The time series of k budgets are even more noisy, so we again apply the 4 h running average and plot smoothed results at the end of the 4 h window.

Figure 6a shows the budget for k between BASE and N165. The actual $\dot{k}$ (from time series of k) becomes negative quickly after seeding starts, then decays towards 0, consistent with the evolution of k for this pair (the dashed blue curve in Fig. 2c). The $\dot{k}$ diagnosed as ${\dot{k}}_L+{\dot{k}}_N$ following Eq. (2) closely matches the actual $\dot{k}$ once we are beyond the first few hours since t_seeding (dashed black line in Fig. 6a). ENTR and SUBS terms drive a more negative k, while the RAD and SURF impose a more positive trend. This qualitative behavior is expected from the l budget in Fig. 5b. Quantitatively, the ENTR term initially dominates but later becomes more comparable in magnitude to RAD, SUBS, and SURF. Encouragingly, the RES term does not bias $\dot{k}$. The contribution of N stays at a steady negative value after initial fluctuation because the N values in BASE and N165 get closer over time, making the negative k even more negative (see Fig. 2b and Eq. 15). In the last few hours of the simulation, ${\dot{k}}_N$ dominates the negative $\dot{k}$, even though the magnitude of ${\dot{k}}_N$ is small compared with the contributions of individual processes via $\dot{L}$. In other words, some of the decrease in k between BASE and N165 in Fig. 2c is due to the fact that BASE continues to gain N but N165 does not (Fig. 2b). The impacts of N evolution are further discussed in Appendix B.

Figure 6 The 4 h smoothed time series for variables related to the k budget, averaged across 22 MCs: (a) for k between BASE and N165 and (b) between N250 and N400.

Figure 6b shows the budget for k between N250 and N400. The qualitative behavior of all terms is similar to those for k between BASE and N165 but with smaller magnitudes.

3.4 Timescales for L and δL

Figure 7a shows the evolution of 22-MC composites of L for BASE, N165, N250, and N400 in the phase space of L–$\dot{L}$. (Note that in this subsection, we use line colors, instead of line types, to indicate aerosol configurations.) The trajectories of all four aerosol configurations start from the circles near (71 g m⁻², 1.8 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) at t_seeding and move towards the ends with no symbols. The circles along the trajectories are 4 h apart. Different from time series shown in Sects. 3.2 and 3.3, these trajectories are based on the 2 min model output and not smoothed. For BASE, the trajectory is roughly linear with L increasing and $\dot{L}$ becoming less positive over time. Therefore, $\dot{L}$ and L can be fit to a linear model:

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\dot{L}=-(L-L_{\mathrm{\infty }})/\mathit{\tau },\end{array}$$

where $\mathrm{1}/\mathit{\tau }$ is the slope and $L_{\mathrm{\infty }}/\mathit{\tau }$ is the intercept. The solution to Eq. (19) is an equation for L exponentially converging towards L_∞ from its initial state L₀ at an initial time t₀:

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }L\left(t\right)=L_{\mathrm{\infty }}+\left(L_{\mathrm{0}}-L_{\mathrm{\infty }}\right)\exp \left(-{\displaystyle \frac{t-t_{\mathrm{0}}}{\mathit{\tau }}}\right)\end{array}$$

(see the solid black curve in Fig. 2a). This evolution manifests in the L–$\dot{L}$ plane as a trajectory approaching (L_∞, 0) and slowing down in that process. Hence, L_∞ can be interpreted as an apparent steady-state L, assuming the linear relation is valid as t approaches infinity; τ can be interpreted as a timescale governing the evolution. These two fitting parameters reasonably describe the evolution within the 24 h simulation period since t_seeding but do not necessarily suggest that the evolution can be characterized by the same equations and/or parameters if the simulations were to last longer. In our case, τ for BASE is between 10 and 15 h (see reference lines in Fig. 7a).

Figure 7 (a) Trajectories starting from t_seeding=12 h, averaged across simulations for BASE, N165, N250, and N400, in the L–$\dot{L}$ plane. All trajectories start around (71 g m⁻², 1.8 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Circles are 4 h apart from t_seeding to 20 h since t_seeding. (b) Same as panel (a) but in the δL–$\mathit{\delta }\dot{L}$ plane, where δ indicates the deviation from BASE. All trajectories start around (0 g m⁻², 0 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). For both panels, gray lines indicate the slopes for reference timescales of 5 h (the steepest one), 10 h, and 15 h (the least steep one).

For three seeded runs, the trajectories deviate from BASE dramatically in a few ways. For N400, the trajectory makes a loop. To be specific, $\dot{L}$ is reduced quickly, becomes significantly negative, and then returns to a positive value. During the same time period, L decreases to around 65 g m⁻². After 8 h since t_seeding, $\dot{L}$ fluctuates around 0.3 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ with L increasing slowly. The trajectories for N165 and N250 behave similarly with smaller loops.

To understand this rather complicated evolution of seeded runs, we decompose their trajectories to the trajectory of BASE and a deviation from BASE, i.e.,

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }L=L_{\text{BASE}}+\mathit{\delta }L,\dot{L}={\dot{L}}_{\text{BASE}}+\mathit{\delta }\dot{L},\end{array}$$

where δ indicates the deviation from BASE. Figure 7b shows the trajectories of three seeded runs in the δL–$\mathit{\delta }\dot{L}$ plane. All three trajectories start from around (0 g m⁻², 0 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) when seeding starts and again move towards ends with no symbols. For N400, $\mathit{\delta }\dot{L}$ quickly reaches its most negative point, after which it oscillates. This oscillation is possibly associated with the redistribution of the seeded aerosol (e.g., Dhandapani et al., 2025) and may be an interesting feature for future investigation. Regardless, it is evident that $\mathit{\delta }\dot{L}$ gradually becomes less negative and approaches zero. As a result, δL becomes more negative but stabilizes. The trajectory approximately follows a line after 4 h since t_seeding. Similarly, we can also define a timescale for this evolution of δL, τ_δL, from the inverse of the $\mathit{\delta }\dot{L}$–δL slope. A rough but reasonable estimate of τ_δL would be 5 h (see reference lines in Fig. 7b).

To understand τ_δL, we compare it with the timescales examined in Jones et al. (2014): a long timescale τ₃ around 3 d, an intermediate timescale τ₂ around 1 d, and a short timescale τ₁ shorter than about 10 h (see their Table 2). The authors attributed τ₃ to the adjustment of z_i and τ₂ to the thermodynamic adjustment in the BL. For τ₁, they resorted to a mechanism called the entrainment–liquid flux (ELF) feedback, which was coined by Bretherton and Blossey (2014) to describe the mechanism for the response in both stratocumulus-topped BLs (STBLs) and cumulus-under-stratocumulus BLs to climate change. For well-mixed STBLs (i.e., conditions explored in our paper), it builds on the strong feedback between w_e, L, and BL turbulence (Zhu et al., 2005): enhanced w_e warms and drys the BL, which reduces L, limits the buoyancy flux, weakens the BL turbulence, and eventually limits w_e. Jones et al. (2014) showed that the signature of this feedback is the relatively steep w_e–ζ_c slope (see their Figs. 2b and 3), which contributes to the short timescale in L evolution via the cloud-base height evolution (see their Eq. 24).

Our τ_δL is clearly comparable to this τ₁. We examine the relation between the 22-MC composites of δw_e and δζ_c for three seeded runs in Fig. 8, where δ indicates the deviation from BASE. All three (δζ_c, δw_e) trajectories start from (0, 0 mm s⁻¹) at t_seeding. Using N400 as an example, δw_e increases and δζ_c decreases until the trajectory reaches its peak after about 2 h. (Note that the solid dots are 2 h apart.) Then, δw_e decreases towards 0 and δζ_c becomes more negative. The negative δw_e–δζ_c slope beyond 2 h is considered “steep” compared with the values of around 15–30 mm s⁻¹ in Jones et al. (2014), suggesting a short timescale via the ELF feedback.

Figure 8 Same as Fig. 7 but in the δζ_c–δw_e plane. All trajectories start from (0, 0 mm s⁻¹) at t_seeding=12 h. Dots are 2 h apart from t_seeding to 12 h since t_seeding.

We further break down δw_e by linearizing Eq. (17):

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }{w}_{\mathrm{e}}=& {\displaystyle }\mathit{\delta }\left({\displaystyle \frac{A}{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}}}\right){\displaystyle \frac{w_{\ast }^{\mathrm{3}}}{z_{\mathrm{i}}}}+\left({\displaystyle \frac{A}{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}}}\right){\displaystyle \frac{\mathit{\delta }\left(w_{\ast }^{\mathrm{3}}\right)}{z_{\mathrm{i}}}}\\ \text{(22)}& {\displaystyle }& {\displaystyle }-\left({\displaystyle \frac{A}{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}}}\right){\displaystyle \frac{w_{\ast }^{\mathrm{3}}}{z_{\mathrm{i}}^{\mathrm{2}}}}\mathit{\delta }{z}_{\mathrm{i}}.\end{array}$$

The three terms on the right-hand side of this equation represent the impacts of responses in $A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}$, $w_{\ast }^{\mathrm{3}}$, and z_i in the seeded runs. This interpretation is valid because all factors other than the variables prefixed with δ are based on the evolution of the BASE and do not correlate with the evolution of δζ_c (not shown). For N400, the $\mathit{\delta }\left(A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}\right)$ term rapidly increases in the first 2 h after seeding starts (dashed green line). During the same time period, the $\mathit{\delta }\left(w_{\ast }^{\mathrm{3}}\right)$ term first increases slightly, then starts to decrease as the enhanced entrainment warming/drying weakens the turbulence. Over the next 2 h, the $\mathit{\delta }\left(A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}\right)$ decreases quickly, likely driven by the quick decrease in N after peaking (Fig. 2b). The trend in $\mathit{\delta }\left(w_{\ast }^{\mathrm{3}}\right)$ is rather weak during this time period. Later, both the $\mathit{\delta }\left(A/\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}\right)$ and $\mathit{\delta }\left(w_{\ast }^{\mathrm{3}}\right)$ terms contribute to the decrease in δw_e. The δz_i term is negative all the time but its magnitude is small. The linearized equation does not fully recover the actual δw_e for N400 (see the gap between the thick green curve and the thin green curve with dots), which is not surprising given the relatively large change in N between BASE and N400. However, the δw_e–δζ_c slope is well captured. N165 and N250 behave qualitatively similar to N400.

To summarize, by decomposing the L in seeded runs into the L_BASE and the deviation from it (δL, which is the L adjustment), we show that δL becomes negative following an exponential convergence. This evolution is governed by a short timescale of 5 h and lasts for 8–12 h. This short timescale is comparable with the one examined in Jones et al. (2014) which the authors attributed to the feedback between w_e, L, and BL turbulence (Zhu et al., 2005). We check the δw_e–δζ_c slope in our simulations to conclude that our fast evolution is likely due to the same mechanism.

3.5 Timescales for k

Since the $\dot{k}$ time series is very noisy, the k timescale is not readily identifiable from trajectories in the k–$\dot{k}$ plane. However, from Eqs. (2), (14), and (15) and results in Fig. 6, one may anticipate that when ${\dot{k}}_L$ dominates, the timescale for k is close to the timescale for l, which is close to τ_δL when the range of L variation discussed here is not too large; later when ${\dot{k}}_N$ becomes important the timescale for k becomes longer because $\dot{k}$ is roughly constant and does not scale with k anymore. To demonstrate this idea we build a two-equation ordinary differential equation (ODE) set for k and δn. One equation is based on Eq. (2),

$$\begin{array}{ll}{\displaystyle }\dot{k}& {\displaystyle }={\displaystyle \frac{\mathit{\delta }\dot{l}}{\mathit{\delta }n}}-k{\displaystyle \frac{\mathit{\delta }\dot{n}}{\mathit{\delta }n}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\dot{l}-{\dot{l}}_{\text{BASE}}}{\mathit{\delta }n}}+{\dot{k}}_N\\ \text{(23)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\delta }n}}\left({\displaystyle \frac{{\dot{L}}_{\text{BASE}}+\mathit{\delta }\dot{L}}{L_{\text{BASE}}+\mathit{\delta }L}}-{\displaystyle \frac{{\dot{L}}_{\text{BASE}}}{L_{\text{BASE}}}}\right)+{\dot{k}}_N,\end{array}$$

and the other equation is rearranged from Eq. (15),

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }\mathit{\delta }\dot{n}=-{\displaystyle \frac{\mathit{\delta }n}{k}}{\dot{k}}_N.\end{array}$$

To close this ODE set, we assume that (1) the evolution of L for BASE is entirely governed by Eqs. (19) and (20) (i.e., with no k or δn in these two equations), (2) δL also exponentially converges (see Fig. 7b), and (3) ${\dot{k}}_N$ is constant (e.g., Fig. 6a). With parameters representing the evolution of k between BASE and N165 in our simulation, the results confirm the reasoning above (Fig. 9).

Figure 9 Evolution of k from the ordinary differential equation (ODE) set for k and n based on Eqs. (2) and (15) with parameters representing the evolution of k between BASE and N165: (a) time series of k and the components from ${\dot{k}}_L$ and ${\dot{k}}_N$ and (b) trajectories of k and its components in the k–$\dot{k}$ plane.

4 Discussion

In this section, we first discuss the timescales associated with L and δL in our simulations and the timescales discovered in previous works. This helps us understand where these timescales fit into a bigger conceptual picture in both the model and the real world. Then we will discuss a few specific issues related to the δL, connecting back to this conceptual picture when necessary.

4.1 Comparison of δL timescales with previous studies

In both Bretherton et al. (2010) and Jones et al. (2014), short timescales emerge in the mixed-layer model (MLM) simulations and LES of stratocumulus as the model states approach their own long-term steady states, which are characterized by the balance between the entrainment velocity and the subsidence velocity in the simulations. (Hereafter we refer to these long-term steady states as the “equilibrium states”, loosely following the terminology in the aforementioned two papers.) In Bretherton et al. (2010), the initial model states are different from the equilibrium states and hence can be interpreted as equilibrium states with perturbations. The short timescale is evident in the first 1 to 2 d of the simulations and decays as the system moves into the slow manifolds governed first by the thermodynamic timescale and later by the inversion adjustment timescale. In Jones et al. (2014), perturbations were introduced into simulations that had reached their equilibrium states, following van Driel and Jonker (2011). After the perturbations, the short timescale dominates for at least 12 h and decays later (e.g., their Fig. 3; also note that several diagnostics used to characterize the short timescale were based on the first 12 h results after perturbation).

In our simulations, the timescale for L in BASE is about 15 h between t_seeding=12 h and the end of the 36 h simulations (Fig. 7a), suggesting that the BASE, after initial spin-up of turbulence, is moving towards the stage dominated by the thermodynamic timescale. So, the states in our BASE simulations are “perturbed” from the equilibrium states for BASE. The seeding introduces another level of perturbation into BASE and the timescales we have examined so far are mostly associated with this perturbation due to seeding. This is different from Bretherton et al. (2010) and Jones et al. (2014), where the reference states are the equilibrium states. Still, after the 30 min seeding ends, δL evolves following a short timescale (Fig. 7b). As in Bretherton et al. (2010) and Jones et al. (2014), the evolutions of our model states are driven by multiple processes, each with its own timescale. However, the fast ELF feedback dominates (Fig. 8) to produce the apparent short timescale. In particular, with the evolution of δN, the k evolution shows varying timescales even before the end of the simulations (Fig. 9).

With this conceptual picture, we can anticipate the timescales found in both our BASE simulations and the deviation of seeded runs from BASE to be transient. Even though it appears that both BASE and seeded runs are approaching some apparent steady states within the 36 h simulations, it takes more than 10 d for them to eventually reach the equilibrium states shown in Bretherton et al. (2010). During this process, these apparent steady states, for example, for L, which are based on the extrapolation of evolution within the simulation period, are likely to drift. From previous works, we expect the evolution of BASE will slow down during this period as the dominant timescale becomes longer due to the change in the primary physical mechanism that drives the evolution. Furthermore, we speculate that the evolution of the deviation from the BASE will also slow down. Hence, for both theoretical understanding and real-world application, it is important to know how long the short timescale dominates. From the composites based on MAIN, δL and hence k evolve according to the short timescale for 8–12 h. For our particular experiment design, since N in some seeded runs (e.g., N400) is still decreasing at the end of the simulation, δL may become less negative as δN values become steady.

4.2 Similarity of δL evolution

In MAIN, even though the magnitude of δL depends on N, the timescale of δL is similar among N165, N250, and N400 (Fig. 7b). As a result, the trajectories in the δL–$\mathit{\delta }\dot{L}$ plane are similar after appropriate scaling. To illustrate this feature, we divide δL by δL^∗, defined as the magnitude of the median δL between 23 and 24 h since t_seeding=12 h; we divide $\mathit{\delta }\dot{L}$ by $\mathit{\delta }{\dot{L}}^{\ast }$, defined as $\mathit{\delta }{L}^{\ast }/\mathrm{5}$ h. With this scaling, a slope of −1 in the $\mathit{\delta }L/\mathit{\delta }{L}^{\ast }$–$\mathit{\delta }\dot{L}/\mathit{\delta }{\dot{L}}^{\ast }$ plane conveniently indicates a timescale of 5 h. The scaled trajectories of N165, N250, and N400 approximately collapse (Fig. 10a). This is also true for similarly scaled δl–$\mathit{\delta }\dot{l}$ (Fig. 10b) as well as scaled δζ_c and scaled δw_e (Fig. 10c). (The scaling factor $\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c}}^{\ast }$ is defined as the magnitude of the median δζ_c between 23 and 24 h since t_seeding; $\mathit{\delta }{w}_{\mathrm{e}}^{\ast }$ equals $\mathrm{20}\mathrm{mm}{\mathrm{s}}^{-\mathrm{1}}\cdot \mathit{\delta }{\mathit{\zeta }}_{\mathrm{c}}^{\ast }$. With this scaling, a slope of 1 in the $\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c}}/\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c}}^{\ast }$–$\mathit{\delta }{w}_{\mathrm{e}}/\mathit{\delta }{w}_{\mathrm{e}}^{\ast }$ plane indicates a δw_e–δζ_c slope of 20 mm s⁻¹.)

Figure 10 (a) Same as Fig. 7b but for scaled (δL, $\mathit{\delta }\dot{L}$). (b) Same as panel (a) but for scaled (δl, $\mathit{\delta }\dot{l}$). (c) Same as Fig. 8 but for scaled (δζ_c, δw_e). (d) Trajectories averaged across 22 MCs in the δζ_c–$\mathit{\delta }{\dot{\mathit{\zeta }}}_{\mathrm{c}}$ plane; all trajectories start around (0, 0 h⁻¹) at t_seeding=12 h. Circles are 4 h apart from t_seeding to 20 h since t_seeding; solid dots are 2 h apart from t_seeding to 12 h since t_seeding.

To further understand this similar τ_δL, we derive a simple model to connect δw_e and δζ_c with simplifications based on the phenomenology in MAIN. We start with the definition of ζ_c in Eq. (13). Assuming z_i does not change with time (termed the “fixed boundary layer depth limit” in Jones et al., 2014),

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{\dot{\mathit{\zeta }}}_{\mathrm{c}}={\displaystyle \frac{{\dot{z}}_{\mathrm{i}}-{\dot{z}}_{\text{cb}}}{z_{\mathrm{i}}}}={\displaystyle \frac{\mathcal{F}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}^{\mathrm{2}}}}.\end{array}$$

Next, we decompose terms in Eq. (25) into a reference state (BASE in our case) and a deviation therefrom due to the seeding (again, prefixed with δ). Considering that only ENTR flux divergences respond strongly to seeding (Fig. 3),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\delta }{\dot{\mathit{\zeta }}}_{\mathrm{c}}={\displaystyle \frac{\mathit{\delta }\mathcal{F}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}^{\mathrm{2}}}}\approx {\displaystyle \frac{\mathit{\delta }{\mathcal{F}}_{\text{ENTR}}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}^{\mathrm{2}}}}=-{\displaystyle \frac{\mathrm{1}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}^{\mathrm{2}}}}\\ \text{(26)}& {\displaystyle }& {\displaystyle }\times \left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}\mathit{\delta }\left(\mathrm{\Delta }{F}_{\langle {\mathit{\theta }}_{\mathrm{l}}\rangle ,\text{ENTR}}\right)+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}\mathit{\delta }\left(\mathrm{\Delta }{F}_{\langle q_{\mathrm{t}}\rangle ,\text{ENTR}}\right)\right).\end{array}$$

Again, due to the weak responses in Δθ_l and Δq_t, we further have

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }{\dot{\mathit{\zeta }}}_{\mathrm{c}}& {\displaystyle }\approx -{\displaystyle \frac{\mathrm{1}}{\langle {\mathit{\rho }}_{\mathrm{0}}{\rangle }_{\text{BL}}{z}_{\mathrm{i}}^{\mathrm{2}}}}\left({\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle q_{\mathrm{t}}\rangle }}\mathrm{\Delta }{q}_{\mathrm{t}}+{\displaystyle \frac{\partial z_{\text{cb}}}{\partial \langle {\mathit{\theta }}_{\mathrm{l}}\rangle }}\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{l}}\right)\mathit{\delta }{w}_{\mathrm{e}}\\ \text{(27)}& {\displaystyle }& {\displaystyle }=c_{\mathrm{1}}\mathit{\delta }{w}_{\mathrm{e}},\end{array}$$

where we denote the prefactor in front of δw_e with c₁ to simplify the notation. Lastly, we approximate the trajectories of (δζ_c, δw_e) after 2 h since t_seeding in Fig. 8 with a linear function,

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }\mathit{\delta }{w}_{\mathrm{e}}=c\left({\mathit{\zeta }}_{\mathrm{c}}-\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c},\mathrm{\infty }}\right)=\mathit{\delta }{w}_{\mathrm{0}}+c\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c}},\end{array}$$

where c refers to the δw_e–δζ_c slope, following the notation in Eq. (22) in Jones et al. (2014), and δζ_c,∞ and δw₀ are parameters for this linear relation that can be interpreted as the apparent steady state of δζ_c and the initial δw_e. Substituting Eq. (28) into Eq. (27) reveals a linear relation between $\mathit{\delta }{\dot{\mathit{\zeta }}}_{\mathrm{c}}$ and δζ_c:

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }\mathit{\delta }{\dot{\mathit{\zeta }}}_{\mathrm{c}}=c_{\mathrm{1}}c\left(\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c}}-\mathit{\delta }{\mathit{\zeta }}_{\mathrm{c},\mathrm{\infty }}\right).\end{array}$$

This is indeed the case with $\mathit{\delta }{\dot{\mathit{\zeta }}}_{\mathrm{c}}$ and δζ_c from the simulations (Fig. 10d). Again, this means δζ_c exponentially converges to its apparent steady-state δζ_c,∞, governed by a timescale $\mathrm{1}/\left(c_{\mathrm{1}}c\right)$.

In Eq. (29), c₁ is mainly controlled by the meteorological conditions (see Eq. 27) in BASE; c describes the phenomenology of the relation between δw_e and δζ_c in our simulations. In a MLM framework, one would be able to show that c is the sum of terms controlled by both meteorology and aerosol (through the entrainment efficiency), e.g., by linearizing w_e in Eq. (15) in Dal Gesso et al. (2014). This exercise is beyond the scope of current work.

Here, we make three points. First, even though the similarity among seeded runs is striking, it is more likely approximate than exact. Based on the conclusions from perturbing the linearized MLM in Jones et al. (2014), the relative changes in the short timescale are expected to directly connect to the relative changes in the entrainment efficiency, which is within 15 % in MAIN (Fig. 4b). With the noisy LES data, this difference may not be evident. If this is the case, the results suggest that the short timescale is mostly determined by the reference state and that the difference between aerosol perturbations within a certain range only slightly modifies it. Second, the results in Jones et al. (2014) suggested that the similarity in c also exists when the equilibrium states in MLM are perturbed (see the similar δw_e–δζ_c slopes in their Fig. 3). For these perturbations, our simplifications (e.g., weak responses in the flux divergences by processes other than entrainment, jumps, and other factors) may not be applicable. It would be interesting to connect the response due to seeding to the response to these perturbations. Both points direct future work to addressing a broader question: how does the short timescale depend on meteorological and aerosol reference states and meteorological and aerosol perturbations? Lastly, as mentioned in Sect. 3.4, δL starts the exponential convergence soon after the seeding stops. After that, the evolution of δL is only controlled by two parameters: the initial $\mathit{\delta }\dot{L}$ (approximated by the most negative $\mathit{\delta }\dot{L}$ along each trajectory in Fig. 7b) and the timescale τ_δL. The similarity suggests that τ_δL is insensitive to the amount of seeded aerosol. Then the evolution of δL (at least for the first 8–12 h) is controlled by a single parameter, the initial $\mathit{\delta }\dot{L}$, which is more negative for greater seeding amount.

4.3 Implications

Proposed MCB efforts involve the seeding of large areas of marine stratocumulus using arrays of aerosol sprayers (Wood, 2021; Feingold et al., 2024). Our results suggest that the L in nocturnal, non-precipitating marine stratocumulus has the highest negative L response per unit amount of N perturbation on a log scale when seeding is applied to cloud fields with low background N, i.e., close to but above the value of N associated with the onset of precipitation. Still, it is unlikely that the negative L adjustment at these values of N will fully offset the brightening due to the Twomey effect during the course of one night (assumed to be 12 h long). With increasing N, the entrainment enhancement and therefore negative L adjustment become less efficient. As a result, stronger N perturbations would likely benefit MCB not only via the Twomey effect but also by reducing the magnitude of the negative adjustment slope. Recently, Prabhakaran et al. (2024) suggested that if non-precipitating marine stratocumuli were to be chosen as targets for MCB, it would be preferable to inject as much aerosol as possible, up until the point that aerosol particle coagulation losses start to dominate. This result is in accord with the results presented in the current work.

The development of the negative δL is a manifestation of the feedback between initially enhanced cloud-top entrainment, L, and BL turbulence soon after the seeding stops. In our simulations, this process could follow a short timescale (around 5 h) for 8–12 h. Our results are robust for STBLs where both the reference and seeded conditions are close to well-mixed and overcast without shear as a primary source of BL turbulence. See Appendices C and D for a few sensitivity tests that support the robustness of our results. It remains to be seen whether additional sources or sinks of BL turbulence, e.g., shear (Kazil et al., 2016) or the presence of a subcloud stable layer (Zhang et al., 2023), and additional timescales introduced by other processes in the real world (e.g., diurnal cycle, drift in SST and FT humidity, spreading of ship tracks) would change the timescale and duration of the δL and k evolution for real clouds.

5 Summary

In this study, we have explored the magnitude and the timescale of liquid water path (L) adjustment to cloud droplet number concentration (N) for nocturnal non-precipitating marine stratocumulus using large eddy simulations (LESs). For each of 22 meteorological conditions (MCs), a 36 h BASE run is performed with relatively low N. Then a set of aerosol perturbation runs is performed by seeding the BASE at t_seeding=12 h and extending the simulations until 36 h. All simulations feature boundary layers (BLs) that are close to well-mixed and overcast, without shear as a primary source of BL turbulence.

The main findings include the following.

• The L adjustment (δL) slope (k) is more negative for simulation pairs with relatively low N and less negative for high N, consistent with Lu and Seinfeld (2005) and Chen et al. (2011). This result agrees with the expectation that the effects of several mechanisms for entrainment enhancement by increasing N all saturate at high N. Overall, a k more negative than −0.4 is unlikely within 24 h since t_seeding.

• The evolution of δL follows a short timescale. To be specific, the δL tendency quickly becomes negative as seeding starts. Soon after seeding ends, it decays towards 0 as δL continues to become negative following approximately an exponential convergence that is governed by a short timescale of around 5 h. The evolution follows this exponential convergence for around 8–12 h.

• Like the short timescale investigated by Jones et al. (2014), our short timescale emerges as a result of the feedback between entrainment, L, and boundary layer (BL) turbulence (e.g., Zhu et al., 2005; Bretherton and Blossey, 2014) driving the δL evolution. To be specific, the seeding enhances the entrainment velocity, which reduces L and weakens the boundary layer (BL) turbulence, leading to the decay of the entrainment velocity deviation in seeded runs from BASE while the responses in entrainment efficiency and integrated buoyancy flux persist. This feedback dominates because other processes, including the radiation, surface fluxes, and subsidence, only respond to the seeding weakly.

• This short timescale is insensitive to the amount of seeding. As a result, the evolution of the deviation of several quantities in seeded runs from BASE is similar after appropriate scaling. The other parameter governing the exponential convergence of δL, namely the initial δL tendency, is sensitive to the amount of seeding.

• The timescale of k evolution is closely related to the timescale of δL and hence also short, while it could also be affected by the timescale of δN.

In Fig. 1, we also observe some correlation between k and meteorological conditions given similar N, like many previous works did. This dependence will be examined in detail in a future study. In the rest of the paper, we have presented results based on multi-MC composite for different aerosol configurations. This is because the time series for individual simulations or individual simulation pairs are noisy, especially when it comes to k and $\dot{k}$. The compositing masks the variation among MCs in terms of both the magnitude and timescale of δL and k evolutions, which should be examined in the future.

Even though LES is advantageous for resolving some small-scale processes, quantifying the effects of individual mechanisms behind the N enhancement of the entrainment using LES takes very careful experiment design (e.g., Igel, 2024) and would be computationally very expensive to do for a wide range of conditions. A companion paper by Hoffmann et al. (2025) addresses this issue in an MLM framework.

Appendix A: Interpretation of k between seeded runs as k due to seeding

We argue that k is more negative for low N and less negative for high N. However, this statement is ambiguous because it takes two cloud states with different N to define a k and it is unclear whether we are talking about the background N or seeded N or some N range that these two states span. Note that in MAIN, all seeded runs are based on the same BASE for each MC. We perform a sensitivity test where we modify the initial N_a in the 22 BASE runs in MAIN so that their N values at 36 h of the simulations are close to 165 cm⁻³. Then we seed these new BASE runs to reach N around 250. We refer to the new BASE as B165 and the new seeded runs as B165N250. The evolution of k between B165 and B165N250 averaged across MCs is quite similar to that between N165 and N250. Based on this results, we can generally claim that k is more negative for lower background N.

Appendix B: Sensitivity to the evolution of N

The evolving N not only directly contributes to k evolution through ${\dot{k}}_N$, but also has a footprint in ${\dot{k}}_L$ via $\dot{L}$ (e.g., Gryspeerdt et al., 2022). In MAIN, N evolves due to surface flux of aerosol, dilution of BL aerosol through mixing with FT air, microphysical processes, and so on. Even though all these processes exist in nature, we perform simulations with approximately constant N to quantify the impacts of the N evolution. We seed all simulations including BASE in MAIN but (1) specify different initial seeded rates and (2) nudge total number concentrations in the BL after initial seeding finishes so that N for all simulations stays about the same as the N at the end of corresponding simulations in MAIN. In these simulations, the k is less negative by about 0.05 for pairs between N165 and simulations with less N. This is consistent with the higher N and hence lower L at the end of the simulations in our constant N simulations than their counterparts in MAIN for aerosol configurations between BASE and N165. The impacts on pairs with larger N are less significant.

Appendix C: Robustness of results

We present three additional sensitivity tests in this section to examine the robustness of our results.

First, our BASE simulations in MAIN are not in equilibrium states and it is unclear what the results would be if we were to let BASE evolve longer before seeding. We therefore redo all seeded runs but delay the seeding time by 12 h so that it is 24 h from the beginning of simulation. The k values between 11 and 12 h after this new seeding time are comparable to the previous k around the same period of time after the original t_seeding (Fig. 1a and c), pointing to the robustness of our results.

Second, in MAIN, there is no nudging or tuning of subsidence to compensate for the radiative cooling in FT. As a result, the FT air could cool down by about 0.5 to 1 K in 36 h simulations. This is not very large compared with our range of initial Δθ_l from 6 to 10 K. Still, we want to make sure that this cooling does not introduce unintended N dependence in simulation results. We therefore rerun the whole ensemble with FT θ_l profiles from 100 m above the z_i and higher nudged to their initial values with a 30 min timescale. This nudging produces simulations with higher L and lower z_i, probably because it suppresses the cloud-top entrainment, but the impacts on k are not evident.

To summarize, the sensitivity to seeding time and nudging of FT θ_l profiles shows the robustness of our results. In both sensitivity tests, τ_δl is shorter than 10 h. This insensitivity is probably because both the BASE and seeded simulations in these sensitivity sets are close to well-mixed and overcast with similar dominant processes for BL turbulence (e.g., with no shear), as in MAIN.

Finally, we test the sensitivity of results to horizontal grid spacing. In MAIN and the three sensitivity sets, we use 200 m horizontal grid spacing, which is relatively coarse. We randomly pick one MC and repeat BASE, N165, N250, and N400 with two finer resolutions, (1) 100 m horizontal grid spacing and (2) 50 m horizontal grid spacing (the latter with a reduced horizontal domain size from 48 to 24 km to save computational resources). For both resolutions, the vertical grid spacing stays unchanged at 10 m. The adjustment slopes produced by both fine-resolution configurations are still more/less negative for lower/higher N but the magnitudes are about 70 %–80 % of those in MAIN. The timescales for δl are too noisy to precisely quantify from a single MC, but they are comparable to that in MAIN.

Appendix D: Impacts of perturbing initial N_a and fixed surface flux

Figure D1 Results from the INT set averaged across 19 MCs: time series of (a) k and (b) L in the context of results from MAIN averaged across the same 19 MCs; trajectories in (c) the $\mathit{\delta }L/\mathit{\delta }{L}^{\ast }$–$\mathit{\delta }\dot{L}/\mathit{\delta }{\dot{L}}^{\ast }$ plane and the $\mathit{\delta }l/\mathit{\delta }{l}^{\ast }$–$\mathit{\delta }\dot{l}/\mathit{\delta }{\dot{l}}^{\ast }$ plane. In panels (c) and (d), all trajectories start from the beginning of simulations; circles are 4 h apart from 4 to 24 h.

Key differences between the design of MAIN in the current work and that of the LES ensemble set used in Glassmeier et al. (2021) (this set hereafter referred to as G21) include (1) interactive surface fluxes in MAIN vs. fixed surface fluxes in G21, (2) SST 0.5 K warmer than the initial surface air temperature in MAIN vs. fixed SST in G21, and (3) seeding at 12 h to achieve regularly spaced N levels in MAIN vs. randomized initial N_a in G21. We simulate one intermediate set (hereafter denoted with INT for “intermediate”) for the 22 MCs used in MAIN to illustrate a few points. This set is configured with fixed surface sensible and latent heat fluxes and a fixed SST that are identical to G21 and the aerosol configuration that is closer to MAIN. For each MC, there is a BASE run using initial N_a that is not greater than BASE in MAIN (L in INT is overall lower than in MAIN and hence can support lower N without precipitating); then there are three perturbation runs where the initial N_a is set up to achieve N values similar to N165, N250, and N400 runs in MAIN. Each simulation lasts for 36 h. One MC in INT has a maximum cloud top hitting the damping layer, while two other MCs have N in BASE higher than the target for N165 at the end of 36 h simulations. These three MCs are removed and we focus on the results from the remaining 19 MCs.

Figure D1a shows the k time series from this INT set. It is more negative for lower N, as in MAIN. In the first 2 h, k between BASE and the seeded run is clearly negative and different between seeded runs, suggesting that N leaves a footprint on the very first few overturnings in the BL, even though the turbulence during that time is not realistic. This is also evident in the L time series in Fig. D1b. However, by 24 h from the beginning of the simulation, k is comparable to those in MAIN at 24 h since t_seeding, suggesting that the less negative k at high N is not a result of compensation for increased entrainment by surface flux, which is consistent with results in Sect. 3 that the response in the surface flux is weak.

The evolution of k is slower in INT than in MAIN, evident in that the k values from INT are less negative than their counterparts from MAIN before they become similar around 24 h. Still, the timescale for k is shorter than 20 h, given that their time series are more curved than reference curves representing exponential convergence of k with a 20 h timescale (Fig. D1a).

Figure D1c shows the (δL, $\mathit{\delta }\dot{L}$) trajectories from INT. After the initial few hours with dramatic fluctuation, δL evolves with a long τ_δL, which arguably becomes shorter after 12 h from the beginning of the simulation. This slow development may be related to the overall weak turbulence in BLs with low L (Fig. D1b) and partial cloudiness (the cloud fraction, f_c, defined as the fraction of domain with cloud optical depth greater than 1, is around 74 % at 2 h from the beginning of the simulation, increases to 92 % at 12 h, and continues to increase towards 36 h; not shown). However, recall that it is the timescale for δl, not δL, that directly relates to the timescale of k. For INT, the timescale for δl is shorter than δL (Fig. D1d), supporting a shorter timescale for k.