Exploring the processes of liquid water path sensitivity to aerosol-cloud interactions using output from a high-resolution large-eddy simulation

Diagnostics from high-resolution Large-Eddy Simulations (LES) are used to investigate aerosol impacts on the liquid water path (LWP) sensitivity in a non-precipitating, single-layer liquid cloud regime. In two LES simulations, the 2013 conditions represent a low aerosol scenario, while the 1985 conditions represent a high aerosol scenario. Joint histograms of cloud droplet number concentration (N_d) and LWP reveal a non-linear relationship, with positive LWP sensitivity (increasing LWP with N_d) at low N_d and negative sensitivity at high N_d. The transition from positive to negative LWP sensitivity occurs at higher N_d values in the 1985 simulation (≈ 300 cm⁻³) compared to the 2013 simulation (≈ 100 cm⁻³), indicating that enhanced aerosol loading shifts the transition point. This shift reflects stronger droplet activation and sustained LWP growth under high cloud condensation nuclei (CCN) conditions. Diagnostics of the cloud dilution indicate that negative LWP sensitivity is linked to enhanced cloud-top entrainment. The temporal evolution of the N_d–LWP relationship confirms increasing dominance of negative sensitivity in the 2013 case, while the 1985 case exhibits weaker LWP depletion. Additionally, aerosol perturbations also influence thermodynamic properties such as the apparent heating/cooling (Q₁) and the moisture sink (Q₂). Specifically, during negative LWP sensitivity phases, stronger cloud-top drying (moisture sinks) is simulated, particularly at high N_d in 2013, consistent with enhanced entrainment/mixing and evaporation-driven cloud dilution. Aerosol perturbations thus modulate both microphysical and thermodynamic processes, producing distinct LWP sensitivity regimes with important implications for understanding aerosol–cloud–climate interactions.

1 Introduction

The aerosol cloud interactions (ACI) and the resulting effective radiative forcing remain a large source of uncertainty when assessing anthropogenic climate change (Forster et al., 2021; Quaas et al., 2020a). The uncertainty in ACI stems from the response of the clouds to the aerosol perturbation (Forster et al., 2020). In liquid clouds, cloud droplets form on an aerosol particle, which can serves as cloud condensation nucleus (Charlson et al., 1992). An increase in atmospheric aerosol leads to an increase in the cloud droplet number concentration (N_d). Twomey (1974) hypothesised that at a constant liquid water path (LWP), an increased aerosol burden leads to clouds with more numerous small droplets, which increase the cloud albedo. In addition, smaller droplets delay the precipitation formation by reducing the collision–coalescence efficiency and increasing the cloud lifetime (Albrecht, 1989). The increase in the response in the N_d also leads to further rapid adjustment of the cloud properties. It includes the alteration of cloud drop size distribution, changes in the LWP, cloud fraction, and dynamic process (Ackerman et al., 2000; Mülmenstädt and Feingold, 2018). Thus, the instantaneous Twomey effect and cloud rapid adjustments contribute to the effective radiative forcing due to ACI (Bellouin et al., 2020; Quaas et al., 2024).

The response of cloud water path, the vertical integral of cloud water, to aerosol perturbation is a key component of cloud adjustments, and yet it is uncertain. This is particularly because of the elusive sign of the LWP adjustment/sensitivity to aerosol perturbations and its regime dependency (Fons et al., 2023; Dipu et al., 2022; Glassmeier et al., 2021; Possner et al., 2020; Gryspeerdt et al., 2019). A positive LWP adjustment is mainly observed in precipitating clouds, in which an increase in the aerosol results in enhanced N_d and smaller droplets, suppressing the precipitation and allowing for an accumulation of LWP (Albrecht, 1989). Thus, the positive LWP adjustment results in thicker and more reflective clouds with a stronger cooling effect. On the other hand, a negative LWP adjustment is associated with cloud droplet evaporation. Aerosol perturbations increase N_d, yielding smaller droplets that evaporate more efficiently during mixing and reduce droplet sedimentation, thereby redistributing liquid water toward the inversion and potentially strengthening cloud-top long-wave radiative cooling. As radiative cooling occurs within a relatively thin layer near the cloud-top, it also promotes the entrainment of warm, dry free-tropospheric air, further enhancing evaporation and potentially reducing LWP - a negative LWP adjustment (Williams and Igel, 2021; Wood, 2012; Bretherton et al., 2007; Ackerman et al., 2004; Wang et al., 2003). Notably, entrainment can reduce LWP even when droplet mass decreases approximately homogeneously, not just due to the preferential evaporation of smaller droplets. Both observational and modelling studies demonstrate a strong offsetting warming effect from negative LWP adjustment (Gryspeerdt et al., 2021; Possner et al., 2020). However, the strength of the net LWP adjustment is modulated by the environmental condition. As a net result of the opposing LWP adjustment mechanisms, their net impact on the large-scale integral remains relatively small or neutral (Zhang et al., 2022). Thus, the bidirectional LWP adjustment/sensitivity, precipitation suppression, and droplet evaporations are difficult to disentangle as these processes coexist in the cloud (Fons et al., 2023).

Recent studies have focused on the sensitivity of LWP to N_d to quantify the impact of aerosols on LWP (Dipu et al., 2022; Gryspeerdt et al., 2019; Bellouin et al., 2020). Using satellite observations, Gryspeerdt et al. (2019) demonstrated that the N_d–LWP relation is non-linear over the global ocean. The LWP adjustment is also regime-dependent (Glassmeier et al., 2021). In marine stratocumulus clouds, the N_d–LWP relationship is non-linear, and the apparent coupling largely reflects co-variability between aerosol loading and meteorological conditions, which fundamentally drives variations in both N_d and LWP (Goren et al., 2025). The LWP adjustment estimated based on satellite observations may be highly uncertain and negatively biased (Arola et al., 2022) because of the retrieval errors and also due to correlated errors in the N_d and LWP retrievals (Quaas et al., 2020b; Gryspeerdt et al., 2019; Grosvenor et al., 2018). In contrast, modelling studies often reported positive LWP adjustments (Quaas et al., 2008; Gryspeerdt et al., 2020). However, high-resolution modelling evidence has also shown negative LWP adjustment, in which the altered N_d leads to enhancement of entrainment mixing, thereby reducing the LWP (Glassmeier et al., 2021; Ackerman et al., 2004). Mülmenstädt et al. (2024a) reported that the latest generation of general circulation models (GCMs) are able to produce negative LWP adjustment besides positive LWP adjustment through precipitation suppression in response to increased N_d. However, the earlier-generation GCMs fail to show negative LWP adjustments in response to anthropogenic aerosols.

Previous studies focused on the N_d–LWP relationship to assist in understanding ACI. Observational studies use natural events where the aerosol perturbation is known and compared with unperturbed cloud regimes (Christensen et al., 2022). Such studies suggest an unchanged LWP or negative LWP adjustment (Malavelle et al., 2017; Toll et al., 2019); however, such cases are limited. Modelling studies, specifically high-resolution Large-eddy simulations (LES), are ideal for improving the understanding of the LWP adjustment by varying aerosol concentration while keeping the other boundary conditions constant. However, the LES simulations are computationally expensive, and the simulations are regime-dependent. Most of the previous LES studies suggest that a positive LWP adjustment is associated with precipitating cloud regimes and a negative LWP adjustment is simulated for non-precipitating cloud regimes (Glassmeier et al., 2021; Ackerman et al., 2004; Hill et al., 2009; Lee et al., 2009; Wang et al., 2011). In this study, we investigate LWP sensitivity in a non-precipitating continental cloud regime using high-resolution LES simulations in numerical weather prediction mode, with initial and boundary conditions from a real weather situation and an interactive land surface (Heinze et al., 2017; Costa-Surós et al., 2020). Aerosol–cloud interaction effects are quantified using control and aerosol-perturbed simulations that differ only in the prescribed cloud condensation nuclei (CCN) fields for droplet activation, while the meteorology is kept identical across the simulations. The following methodology section describes the model setup and aerosol perturbation. Using the same LES simulations, Dipu et al. (2022) demonstrated that the N_d–LWP sensitivity is bidirectional. Using the perturbed aerosol simulation of the same cloud regime, here we investigate the impact of aerosol on bidirectional LWP sensitivity. This allows us to infer the degree to which the N_d–LWP relationship represents a causal influence of N_d on LWP. We examine the impact of aerosol on positive and negative LWP sensitivity and investigate the microphysical and thermodynamic processes controlling the sign and magnitude of LWP sensitivity.

2 Data and Methodology

LES using the ICOsahedral Nonhydrostatic (ICON) model (Dipankar et al., 2015; Zängl et al., 2015) are analysed in this study. The atmospheric model ICON has been configured to a large-eddy simulation framework (Dipankar et al., 2015), and has been evaluated against standard LES models and multiple observations (Heinze et al., 2017). The high resolution, ICON-LES simulation has been performed as part of the High Definition Clouds and Precipitation for advancing Climate Prediction (HD(CP)²) project. The simulation ran over a large domain (over Germany) in a weather prediction mode, which uses realistic boundary conditions from the operational COSMOS-DE (Consortium for Small Scale Modelling, Baldauf et al., 2011), with a fully interactive land surface (Costa-Surós et al., 2020). The model is configured with a horizontal resolution of 156 m and 150 vertical levels with a model top at 21 km. A sub-grid scale turbulence scheme based on the classical Smagorinsky scheme is used in the model, which also accounts for thermal stratification (Lilly, 1962). The model uses a two-moment liquid and ice-phase bulk microphysics scheme (Seifert and Beheng, 2006). In the two-moment microphysical scheme the grid-scale cloud droplet nucleation rate is estimated as a function of CCN concentration, vertical velocity, and supersaturation (Eq. 7 of Seifert and Beheng, 2006). Following Kogan and Martin (1994), the two-moment scheme applies the standard saturation adjustment technique to treat condensational growth of cloud droplets. The CCN concentrations in the model are prescribed as a spatially and temporally varying distribution. The control simulation uses CCN concentrations as estimated for 2 May 2013 (Costa-Surós et al., 2020), and for the perturbed simulation, CCN concentrations valid for the year approximately 1985 were selected, in which the pollution level in Europe was at its peak (Smith et al., 2011). The 2013 CCN concentrations are generated from 2013 emissions using a regional coupled model system (Wolke et al., 2004, 2012). The 1985 CCN concentrations are obtained by scaling the 2013 CCN concentrations with species-dependent factors derived from emission ratios following Genz et al. (2020). A detailed description is provided in Costa-Surós et al. (2020). The simulations were performed over Germany for selected dates, of which the date 2 May 2013 is considered in the study based on the evaluation results from Heinze et al. (2017). The 2 May 2013 has been one of the extensive measurement campaigns for the HD(CP)² Observational Prototype Experiment (HOPE, Löhnert et al., 2015; Madhavan et al., 2016) and the evaluation results from Heinze et al. (2017) suggest the presence of a wide range of cloud regimes, compared to other HD(CP)² simulations. A detailed description of the ICON-LES model and HD(CP)² simulations can be obtained from Dipankar et al. (2015), Heinze et al. (2017), and Costa-Surós et al. (2020).

Here, we have used the coarse-gridded data with a resolution of 1.2 km (grid size of 589 × 637), a standard reduced-volume product that has been used in previous studies and evaluations (Costa-Surós et al., 2020; Dipu et al., 2022). The actual ICON-LES simulation was performed with 156 m horizontal resolution, though. Our results rely on regime-conditioned, cloud-only statistics that mitigate grey-zone smoothing. While coarse-gridding may influence quantitative values, the qualitative sensitivities remain robust. Five-minute instantaneous model output from 1000 to 2000 h is considered for the study. The analysis is restricted to single-layered liquid clouds by excluding the clouds with a cloud-top temperature below 273 K. Cloud-top is defined as the uppermost model level with liquid cloud water present (liquid water content, LWC $>\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$ kg kg⁻¹). We additionally constrained N_d>2 cm⁻³ and restricted the analysis to overcast and optically detectable clouds (cloud fraction = 1 at 1.2 km and cloud optical thicknesses greater than 2) to minimise cloud-edge contamination. In this study, N_d is primarily intended as a consistent metric for comparing aerosol perturbations and linking to cloud-top processes (radiative cooling and entrainment) that control LWP adjustments. Therefore, the cloud-top N_d is not interpreted as a volume-mean value but as an indicator of the droplet population near the region most relevant for the LWP adjustment mechanisms. In the specific ICON-LES, entrainment is not parameterised; it arises from resolved advection and subgrid turbulent mixing. In the selected case, the cloud-top height is at ≈ 1700–2000 m, with an effective vertical resolution of ≈ 50–80 m. Because of the high model resolution and high frequency of model output, the N_d is divided into logarithmic bin sizes of 1000 numbers. The corresponding bin mean cloud microphysical properties are used for the analysis. Additionally, to quantify the grid scale impact of aerosol, cloud properties at the same grid points for both simulations are considered, assuming that the initialization of the cloud fields leads to approximately in the same location in both simulations.

The apparent heating (Q₁) and the moisture sink (Q₂) at the cloud-top are calculated by following Yanai et al. (1973). Q₁ represents the apparent heating/cooling of the atmospheric layer due to various processes such as radiation, condensation, and convection. The corresponding equation is given by,

$$\begin{array}{}\text{(1)}& Q_{\mathrm{1}}=c_{\mathrm{p}}{\displaystyle \frac{\partial T}{\partial t}}-c_{\mathrm{p}}(\mathit{\omega }\mathit{\sigma }-V\cdot \mathrm{\nabla }T)\end{array}$$

where, $\mathit{\sigma }=\left(\frac{RT}{c_{\mathrm{p}}P}\right)-\left(\frac{\partial T}{\partial p}\right)$, the static stability, p is the pressure, V is the horizontal velocity vector, ∇ horizontal gradient operator, c_p is the specific heat of dry air at constant pressure, T is the temperature, ω is the vertical p velocity, and t is the time. Additionally, Q₂ represents the drying/moistening due to condensation or evaporation and moisture flux convergence, and it is represented as,

$$\begin{array}{}\text{(2)}& Q_{\mathrm{2}}=-L{\displaystyle \frac{\partial q}{\partial t}}-LV\cdot \mathrm{\nabla }q-L\mathit{\omega }{\displaystyle \frac{\partial q}{\partial p}}\end{array}$$

where q is the specific humidity and L is the latent heat of condensation. The adiabatic fraction (f_ad) is a measure of cloud dilution, which is primarily due to entrainment, turbulent mixing and evaporation, and is defined as the ratio of LWP to adiabatic LWP (LWP_ad), and is expressed as

$$\begin{array}{}\text{(3)}& f_{\mathrm{ad}}={\displaystyle \frac{\mathrm{LWP}}{{\mathrm{LWP}}_{\mathrm{ad}}}}\end{array}$$

Where, LWP is the liquid water path and LWP_ad is the adiabatic liquid water path (see Appendix A). The joint histograms analysed in this study are constructed as conditional probabilities (CP [%]) following Gryspeerdt et al. (2016) and are defined as the probability of finding a certain LWP given that a certain N_d has been observed (CP = [P (LWP ∣ N_d) × 100 %]). For joint histogram analysis, the variables are binned with a bin size (number of bins) of 1000. In the following analysis, N_d-bin mean variables are used, which means $\stackrel{\mathrm{‾}}{\text{variable}}$ at certain N_d bins (P (variable ∣ N_d) ). Because the joint histograms use the conditional probability, they provide a regime-conditioned, distribution-based statistics that is less sensitive to the exact spatial co-location of individual clouds. In a large LES domain, environmental heterogeneity and aerosol–meteorology co-variability can still influence the apparent N_d–LWP relationship. However, our experimental design compares a control and an aerosol-perturbed simulation with identical meteorological forcing, differing only in the imposed CCN fields. Consistent with the HD(CP)² approach, Costa-Surós et al. (2020) showed that aerosol signals can be quantified in such large-domain LES using domain-wide, regime-conditioned statistics rather than pointwise cloud matching. Accordingly, we use regime-conditioned statistics and distribution-based metrics to mitigate meteorological confounding while retaining the aerosol signal/contrast imposed by the CCN perturbation. Even with the identical large-scale forcing, the vertical velocity w would respond if the aerosol perturbation changes cloud-top cooling/evaporation and precipitation sufficiently. So, in this particular simulation, some dynamical differences (in vertical velocity) are expected in at least parts of the domain/time.

3 Results

The aerosol perturbation in the 1985 simulations results in a significant increase in the number of cloud droplets, as demonstrated by the shift in the N_d probability density function (PDF) distribution (Fig. 1a) towards higher N_d. In 2013, the N_d is primarily distributed between 2 and 700 cm⁻³. In contrast, the 1985 simulation shows a broader distribution that extends up to 1000 cm⁻³, suggesting greater activation of N_d due to aerosol perturbation. Furthermore, this perturbation results in a 120 % increase in the mean N_d compared to the 2013 simulation. However, the LWP distribution shows relatively small shifts towards higher LWP in the 1985 simulation compared to 2013 (Fig. 1b). The relatively small shift in 1985 indicates that aerosol perturbation has less impact on the LWP, contributing only to a 5 % increase in the mean LWP when compared to the 2013 simulation. The contrasting response of the N_d and LWP to the aerosol perturbation suggests that higher droplet activation alone does not directly translate into proportional increases in bulk water content, highlighting the importance of compensating microphysical and dynamical processes (e.g., entrainment/mixing and cloud dilution) that can offset LWP increases.

Figure 1 Comparison of (a) N_d (cm⁻³) and (b) LWP (gm⁻²) probability density function (PDF) for the 2013 and 1985 simulations. The green line denotes the 2013 simulation using present-day (2013) CCN concentrations, while the 1985 simulation applies CCN concentrations representative of peak aerosol loading over Europe around 1985.

To disentangle these processes, we have extended the analysis to the joint histogram of LWP and N_d. Figure 2 shows the joint histogram between LWP and N_d for the 1985 and 2013 simulations. In both cases, the maximum CP depicts a spread for LWP values ranging from 10 to 500 gm⁻², particularly at low N_d (< 100 cm⁻³). At the middle N_d values (between 100 and 500 cm⁻³), the CP narrows, with the highest values occurring in this range for both simulations. At higher N_d, the CP spread increases further, especially in 2013, though CP values above 700 cm⁻³ remain low. This low CP shows that there are limited cloud regimes with high N_d in the 2013 simulation. It is also evident in the N_d distribution (PDF) shown in Fig. 1a.

Figure 2 The N_d-LWP joint histogram for the (a) 2013 and (b) 1985 simulations. The thick black line in each plot shows the smoothed mean LWP ($\stackrel{\mathrm{‾}}{\text{LWP}}$) at certain N_d bins (P (LWP ∣ N_d)). CP(%) is condition probability: the probability of finding a certain LWP given certain N_d. The green line denotes the 2013 simulation using present-day (2013) CCN concentrations, while the 1985 simulation applies CCN concentrations representative of peak aerosol loading over Europe around 1985.

The N_d-bin mean LWP ($\stackrel{\mathrm{‾}}{\text{LWP}}$) in the joint histogram implies the tendency of the N_d-LWP relation. For both simulations, the $\stackrel{\mathrm{‾}}{\text{LWP}}$ shows a non-linear relationship consistent with prior studies (Gryspeerdt et al., 2019; Dipu et al., 2022), but the magnitude and the N_d threshold of the transition differ due to differences in meteorological regime, sampling/aggregation scale, and aerosol loading. This non-linear relation implies a positive LWP sensitivity ($\stackrel{\mathrm{‾}}{\text{LWP}}$ increasing with N_d) for lower N_d values and a negative LWP sensitivity ($\stackrel{\mathrm{‾}}{\text{LWP}}$ decreases with increasing N_d) for higher N_d values. The transition point where the positive LWP sensitivity shifts to negative varies with the aerosol perturbation. In the 2013 simulation, the transition from positive to negative LWP sensitivity is simulated around N_d ≈ 100 cm⁻³ and in the 1985 simulation, it shifts to N_d ≈ 300 cm⁻³. This shift toward higher N_d reflects higher CCN concentrations in 1985, which enhanced droplet activation during cloud formation, increasing N_d and delaying LWP depletion. Since the two simulations have identical meteorological forcing and initial conditions, and differ only in CCN, the shift in the transition point is attributed to the aerosol perturbation rather than thermodynamic warming. The effect of aerosol perturbation is also evident in the R_eff. In the 2013 simulation, $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}}$ peaks at low N_d and then declines sharply beyond ≈ 100 cm⁻³ indicating stronger cloud depletion. Conversely, in the 1985 simulation, the $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}}$ decreases at larger N_d ≈ 300 cm⁻³ indicating thicker, less diluted clouds (figure not shown). The adiabatic fraction (f_ad), defined as the ratio of LWP to the adiabatic LWP (LWP_ad), is a measure of dilution/subadiabaticity, which is primarily due to entrainment, turbulent mixing, and evaporation. Values of f_ad ≈ 1 indicate near-adiabatic, weakly diluted clouds (typically cloud cores), whereas f_ad ≪ 1 denotes strongly diluted, subadiabatic conditions, most commonly near cloud-top and cloud edges. In general, f_ad can increase with N_d at low to moderate N_d when enhanced droplet number delays collision–coalescence and reduces precipitation loss, allowing liquid water to accumulate. At higher N_d, the droplets are smaller and tend to evaporate more efficiently during entrainment/mixing, especially near the cloud-top, leading to cloud dilution and reducing f_ad. Figure 3 shows the relation between bin mean $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$ and the $\stackrel{\mathrm{‾}}{f_{\mathrm{ad}}}$. In both simulations, the lower $\stackrel{\mathrm{‾}}{f_{\mathrm{ad}}}$ represent subadiabatic cloud regime. For positive LWP sensitivity, $\stackrel{\mathrm{‾}}{f_{\mathrm{ad}}}$ increases slightly with increasing $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$ in both simulations. At higher $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$, in the 1985 simulation, $\stackrel{\mathrm{‾}}{f_{\mathrm{ad}}}$ continues to increase with $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$, particularly during positive LWP sensitivity, indicating a relatively weaker dilution (i.e., more adiabatic). However, in the 2013 simulation, the $\stackrel{\mathrm{‾}}{f_{\mathrm{ad}}}$ shows a decrease at higher $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$, suggesting that a more subadiabatic cloud regime implies cloud dilution driven by entrainment/mixing. The analysis further suggests that in the 1985 simulation, the aerosol perturbation leads to thicker, higher LWP clouds at high N_d with weaker cloud-top dilution, whereas the 2013 simulation exhibits thinner clouds at high N_d, indicative of more effective entrainment-driven dilution.

Figure 3 The N_d-bin mean adiabatic fraction (f_ad, %) for the 2013 and 1985 simulations. The points represent the $\stackrel{\mathrm{‾}}{f_{\mathrm{ad}}}$ at certain P (f_ad ∣ N_d ). The diamond and the circle shape denote the corresponding mean LWP ($\stackrel{\mathrm{‾}}{\mathrm{LWP}}$) for 2013 and 1985 simulations, respectively.

Furthermore, the temporal evolution of the cloud field in both simulations illustrates time-dependent LWP sensitivity during the simulations. Figure 4 illustrates the temporal evolution of the LWP sensitivity (time evolution of the joint histogram between $\stackrel{\mathrm{‾}}{\mathrm{LWP}}$ and N_d) for both simulations. Initially, both simulations exhibit positive LWP sensitivity because CCN acts to elevate N_d, reduce droplet size, and suppress or delay warm-rain formation, allowing LWP to accumulate (Albrecht, 1989). Both simulations use identical initial conditions and forcing and differ only in the prescribed CCN, so this positive LWP sensitivity reflects rapid microphysical adjustment rather than differences in initialisation. Over time, however, the relationship becomes more non-linear, with an increasing contribution from negative sensitivity. In the 1985 simulation, the positive LWP sensitivity dominates due to the high CCN concentration, which accounts for the activation of numerous smaller droplets. Additionally, the transition of positive to negative LWP sensitivity shifts toward a much higher N_d value in the 1985 simulation over time. Glassmeier et al. (2021) reported that the LWP sensitivity becomes increasingly negative over time in stratocumulus clouds. Likewise, the LWP sensitivity becomes more negative in both simulations over time. The numerical value of the negative LWP sensitivity is derived as the slope of the linear regression through the $\stackrel{\mathrm{‾}}{\mathrm{LWP}}$ in the specific N_d-bins, depicted in Fig. 5. Notably, the 2013 simulation exhibits a steeper negative slope than the 1985 simulation, indicating a more rapid LWP depletion. Although the transition to negative LWP sensitivity shifts to higher N_d in 1985, the magnitude of the negative slope is smaller than in 2013. The steeper negative slope in 2013 suggests that the cloud field enters drying/entrainment-driven depletion more readily, leading to more rapid LWP depletion. In contrast, the 1985 perturbation primarily shifts the transition to higher N_d and delays the onset of the depletion regime. Importantly, a higher aerosol loading does not necessarily imply a more negative N_d–LWP slope. The strength of the negative LWP sensitivity is state-dependent and depends on whether the clouds are in a dilution-dominated (entrainment–evaporation) state.

Figure 4 Temporal evolution of the mean LWP ($\stackrel{\mathrm{‾}}{\text{LWP}}$) at certain N_d bins (P (LWP ∣ N_d) ) for (a) 2013 and (b) 1985 simulations. Each line indicates the N_d-LWP relationship at every model time step (5 min interval), and the colour gradient indicates the temporal evolution of the N_d-LWP relationship (from 10:00 UTC to 20:00 UTC).

Figure 5 The magnitude of negative LWP (g m⁻²) adjustment, as calculated by the N_d-LWP slope over time. The green diamond shape denote 2013 simulation, and the red circles denote the 1985 simulation. The respective dotted line indicates the linear regression.

From the temporal evolution of the LWP sensitivity, the critical $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$ and $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}}$ values are estimated. The critical values are specified only for time steps at which the domain-wide $\stackrel{\mathrm{‾}}{\mathrm{LWP}}$ exhibits a non-linear dependence on N_d, with a clear maximum (i.e., a transition from positive to negative LWP sensitivity). From the time the temporal evolution of the LWP sensitivity (Fig. 4), if a non-linear relationship is present, the mean N_d $/$ R_eff of the bin at which LWP reaches its maximum (sign change in sensitivity) is taken as the critical value for that time step. If no clear non-linear relationship exists, no critical values are assigned. They are derived across the same spatial domain and time interval. Figure 6 illustrates the critical $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$ and the corresponding $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}}$ at the transition between positive and negative LWP sensitivity. In the 2013 simulation, the critical $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$ is mostly distributed between 48 and 490 cm⁻³, with maximum density centred around 100 to 200 cm⁻³, whereas in 1985, it ranges from 140 to 711 cm⁻³, with the maximum density centred around 400 cm⁻³ in the 1985 simulation. This confirms the rightward shift in the LWP transition point under high-aerosol conditions. At this critical $\stackrel{\mathrm{‾}}{N_{\mathrm{d}}}$, $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}}$ shows a relatively widespread distribution, between 4.7 and 8 µm , in the 2013 simulation. Meanwhile, the 1985 simulation shows a relatively small spread in the $\stackrel{\mathrm{‾}}{R_{\text{eff}}}$ distribution, between 4.2 and 6.5 µm further reflecting the more variable cloud dilution in the low-aerosol case.

Figure 6 Violin plots for critical (a) N_d (cm⁻³) and (b) R_eff (µm). The critical indicates the N_d $/$ R_eff at which the LWP adjustment becomes negative over time. The green colour denote 2013 simulation, and the red colour denotes the 1985 simulation. On each side, the grey line indicates the distribution shape of the data (PDF). The white dot on the violin plot represents the median, the back bar in the centre represents the interquartile range (first and third quartile), and the lower and upper parts of the violin plot represent the lower/upper adjacent values.

Further insights into aerosol-induced cloud changes are revealed by analysing thermodynamic diagnostics, particularly the apparent heating (Q₁) and the moisture sink (Q₂) at the cloud-top. The Q₁ in a cloud layer is associated with various processes, local temperature changes, advection, and vertical motion. A substantial part of apparent heating/cooling is also associated with cloud microphysical processes in clouds. In liquid clouds, condensation of water vapour contributes to microphysical heating, and cloud droplet evaporation contributes to cooling (Keshtgar et al., 2023). Figure 7a shows the relationship between $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ (N_d-bin mean) and N_d at the cloud-top. For lower N_d values, the $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ is negative, implying apparent cooling in both simulations. A negative $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ is mainly associated with cloud droplet evaporation, rising motion, and cold air advection in a cloudy layer. Specifically, in the selected case, at lower N_d, the $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ is negative, though it becomes less negative as N_d increases. The negative $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ observed at lower N_d may be related to processes such as droplet evaporation and/or rising motion. As the cloud develops, collision–coalescence shifts the droplet spectrum toward larger $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}}$ and promotes precipitation processing. As a result, the relative contribution of evaporative cooling associated with mixing/entrainment decreases, and $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ becomes less negative with increasing N_d. This is consistent with concurrent increases in LWC and decreases in specific humidity at lower N_d (Fig. A2 in the Appendix). In the 2013 simulation, as N_d increases (> 200 cm⁻³), $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ becomes positive, indicating apparent heating. However, this apparent heating is only simulated at higher N_d (> 800 cm⁻³) in the 1985 simulation. A positive $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ is associated with processes such as condensation and latent heat release, sinking motion, and entrainment/mixing of warm, dry air. In the 2013 simulation, however, the specific humidity increases, and the water content (LWC/LWP) decreases as N_d > 200 cm⁻³. Therefore, the only possible mechanism left to explain the positive $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ is the entrainment/mixing of warm, dry air and the resulting cloud droplet evaporation, which is in agreement with the dilution of clouds with higher N_d (Fig. 3). Similarly, in the 1985 simulation, the negative $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ further decreases (less negative) for higher N_d (> 300 cm⁻³) and a positive $\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$ is simulated only at higher N_d (> 800 cm⁻³), suggesting that aerosol perturbation shifts the onset of entrainment-induced cloud depletion to higher N_d.

Figure 7 The N_d-bin mean (a) Q₁ and (b) Q₂ for the 2013 and 1985 simulations. The solid lines represent the smoothed mean of the mean Q₁ ($\stackrel{\mathrm{‾}}{Q_{\mathrm{1}}}$) and mean Q₂ ($\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$) at certain N_d bins (same as Fig. 1). The shaded region represents the rolling standard deviation of the respective N_d-bin mean values. The green line denotes the 2013 simulation using present-day (2013) CCN concentrations, while the 1985 simulation applies CCN concentrations representative of peak aerosol loading over Europe around 1985.

The apparent moisture sink, Q₂ (N_d-bin mean, $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$), for the two simulations is depicted in Fig. 7b, which also supports the above interpretation. Generally, a positive Q₂ indicates moisture removal through condensation or dry air advection, while a negative Q₂ implies moisture addition through evaporation or moist air advection. In the 2013 simulation, $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$ is negative for the lower N_d (< 100 cm⁻³), indicating the dominant influence of moist air advection, along with cloud dilution. Simultaneously, the $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$ increases as N_d increases for lower N_d, suggesting a reduction in cloud dilution, which correlates with an increase in LWC and a decrease in specific humidity (Fig. A1 in the Appendix). In contrast, for higher (N_d >100 cm⁻³), in the 2013 simulation, $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$ tends to be positive, indicating moisture removal, consistent with entrainment-driven evaporation. The 1985 simulation shows predominantly positive $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$, except for very low N_d (> 10 cm⁻³). The positive $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$ for N_d values less than 300 cm⁻³ is driven by condensation. At higher N_d (> 300 cm⁻³), however, the same positive $\stackrel{\mathrm{‾}}{Q_{\mathrm{2}}}$ results from cloud water removal through dry air entrainment. This is also consistent with a concurrent increase in specific humidity and a decrease in LWC, indicating cloud dilution associated with entrainment.

4 Discussion

The ambiguity in the LWP adjustment/sensitivity due to aerosols varies with individual clouds (Smalley et al., 2024; Christensen et al., 2020), which adds uncertainty to effective radiative forcing due to the aerosol-cloud interactions (Mülmenstädt and Feingold, 2018). Recent studies utilise the sensitivity of N_d to LWP to improve the understanding of the aerosol-cloud interaction using modelling and observations (Gryspeerdt et al., 2019; Dipu et al., 2022; Mülmenstädt et al., 2024b), in which N_d represents an indirect proxy for aerosols. In a specific cloud regime, a positive LWP sensitivity implies a systematic increase of LWP with increasing aerosols, and a negative LWP sensitivity indicates cloud depletion. Here, we have investigated the significance of aerosol perturbation to the LWP sensitivity using the ICON-LES model.

In the selected cloud regime, the 1985 simulation, representing a high-aerosol scenario, consistently showed systematically higher N_d values under identical meteorological conditions compared to the 2013 simulation (low-aerosol scenario). The high CCN in the 1985 simulation led to more numerous and smaller droplets, sustaining positive LWP sensitivity compared to 2013. As a result, the transition point from positive to negative to positive LWP sensitivity has shifted to a higher N_d (>300 cm⁻³), compared to that (>100 cm⁻³) in 2013. This shift in the negative LWP sensitivity is also visible in the time evolution of the LWP sensitivity in the respective simulation. Despite this shift, both simulations exhibited negative LWP sensitivity at higher N_d, linked to cloud dilution due to entrainment. Analysis of the thermodynamic diagnostics further reveals that in both simulations, cloud-top Q₁ (apparent heating) and Q₂ (moisture sink) became increasingly positive with N_d, indicating warm, dry air entrainment and associated evaporation. The positive trend is more pronounced in 2013, consistent with stronger dilution. Meanwhile, in 1985, Q₁ and Q₂ remained less positive, suggesting a more buffered response to entrainment due to sustained cloud development under high aerosol loading. Additional diagnostics of the temperature tendency term in Q₁ (Eq. 1) at cloud-top reveal that it becomes increasingly negative as N_d increases, particularly for negative LWP sensitivity (Fig. A2a). The positive advection terms (the sum of horizontal and vertical advection) in Q₁ (Eq. 1) for low N_d indicating net warming due to advection, which later transitioned to neutral or negative values at higher N_d, indicating reduced warming and as entrainment increases (Fig. A2c). The high CCN in the 1985 simulations resulted in a less negative temperature tendency (less cooling) and advection-induced cooling compared to the 2013 simulation. Similarly, specific humidity tendencies term in Q₂ (Eq. 2) also showed increased moisture loss at higher N_d, particularly in 2013, further confirming the role of entrainment-driven drying. In the 2013 case, the net moisture advection term in Q₂ (Eq. 2) is initially positive at lower N_d, and it becomes negative at higher N_d, consistent with entrainment-induced drying (Fig. A2b). Conversely, in the 1985 simulation, net moisture advection remains persistently negative and intensifies slightly with increasing N_d, suggesting a weak drying tendency under high aerosol conditions (Fig. A2d).

The analysis shows that the aerosol perturbation induces a clear microphysical and thermodynamic response and modulates the regime-conditioned N_d–LWP sensitivity, while the change in the bulk LWP distribution remains modest. In the 2013 simulation, Q₁ becomes less negative and positive for negative LWP sensitivity (N_d > 100 cm⁻³). Here, entrainment refers to turbulent mixing of relatively warm, dry free-tropospheric air across the inversion into the cloud-top layer, which promotes subsaturation and evaporation. Since Q₁ is an apparent heating residual that includes advection and vertical motion, its sign reflects the net balance of warming and evaporative cooling rather than evaporation alone. Moreover, this evaporation, along with a decrease in local specific humidity, contributes to a moisture sink (Q₂ > 0), particularly for the negative LWP sensitivity. However, in the 1985 simulation, aerosol perturbations led to a less positive tendency in Q₂ for negative LWP sensitivity compared to prior simulations. This can be further explained by the weaker temperature tendencies, reduced local specific humidity, and decreased moisture sink in the 1985 simulation for negative LWP sensitivity. Furthermore, in the 1985 simulation, condensation processes dominate (positive Q₂), which helps maintain a positive LWP sensitivity even at higher N_d. Nonetheless, as N_d increases further, Q₂ also increases, which may be due to enhanced droplet evaporation linked to warm air entrainment. It should be emphasised, however, that evaporation and advection are not the only processes associated with negative LWP sensitivity. While they appear to be the dominant mechanisms in the above simulations, additional processes such as cloud-top radiative cooling, droplet sedimentation, and turbulence entrainment feedbacks may also play important roles in driving cloud depletion at high N_d. We also note that numerical diffusion can affect the magnitudes of the diagnosed apparent heating and moisture sink near cloud edges, which is a limitation of the present analysis and should be considered when interpreting the results. Future studies should aim to disentangle and quantify the relative contributions of these pathways, for example, by combining targeted LES experiments with process-level diagnostics and Lagrangian cloud tracking.

5 Conclusions

This study uses the ICON-LES model as part of the HD(CP)² project to investigate the effect of aerosols on LWP sensitivity. Simulations were conducted over Germany on 2 May 2013, with high (1985 CCN condition) and low aerosol (2013 CCN condition) scenarios. The joint histogram analysis reveals a non-linear relationship between N_d and LWP in both simulations, which is consistent with previous studies. The non-linear relationship implies for low N_d values, LWP increases with N_d (positive LWP sensitivity), while at higher N_d, LWP decreases (negative LWP sensitivity). The transition from positive to negative LWP sensitivity occurs at a lower N_d (≈ 100 cm⁻³) in the 2013 simulation and shifts to higher N_d (≈ 300 cm⁻³) under the 1985 aerosol scenario. This indicates that increased aerosol concentration leads to sustained droplet activation, thereby shifting the cloud depletion to higher N_d. The 1985 simulation exhibits more persistent positive LWP sensitivity, associated with enhanced droplet activation and thicker clouds. In contrast, the 2013 simulation reveals a greater degree of cloud dilution, as indicated by a more pronounced decrease in the R_eff and f_ad at higher N_d. Temporal analysis of the N_d-LWP joint histogram further illustrates a non-linear relationship with negative LWP sensitivity becomes dominant over time in 2013. However, in the 1985 simulation, the positive LWP sensitivity is dominant, with weaker negative LWP sensitivity observed over time.

Furthermore, thermodynamic features such as cloud-top apparent heating (Q₁) and moisture sink (Q₂) also reveal a significant impact of aerosol perturbation. Our analysis suggests that negative Q₁ dominates at low N_d, due to droplet evaporation and/or rising motion. In contrast, the apparent heating (positive Q₁) observed at higher N_d is attributed to cloud dilution and warm air advection. In the 2013 simulation, we found apparent heating or positive Q₁ for N_d > 200 cm⁻³. In contrast, the 1985 simulation showed positive Q₁ only at a much higher N_d value (> 800 cm⁻³). Thus, the aerosol perturbation results in sustained negative Q₁ for higher N_d, with a weaker positive Q₁ or cloud dilution through entrainment. Similarly, negative Q₂ or moisture gain is simulated at low N_d despite the negative Q₁, indicating dominant cloud growth. While positive Q₂ or moisture sink is simulated at higher N_d, indicating a drying effect through warm air entrainment. Our analysis suggests that high CCN concentration in the 1985 simulation exhibits greater moisture retention (negative Q₂), supporting sustained cloud growth and positive LWP sensitivity to higher N_d. In contrast, the moisture sink (positive Q₂) is observed in relatively higher N_d (> 800 cm⁻³) in the 1985 simulation.

Both simulations reinforce the hypotheses that negative LWP sensitivity at high N_d are closely associated with entrainment-driven cloud dilution, evidenced by increased Q₁ (warming), increased Q₂ (moisture loss), reduced LWC, and droplet evaporation. However, the threshold N_d for dilution shifts to higher N_d in the 1985 simulation, indicating the enhanced effect of aerosol perturbation on mitigating cloud depletion. The response in the N_d–LWP relationship under aerosol perturbations implies a modified cloud radiative response, with sustained positive LWP sensitivity enhancing cloud albedo towards high N_d with weaker cloud depletion through entrainment at higher N_d. These dynamics are critical for quantifying the effective radiative forcing of aerosol-cloud interactions in convective cloud regimes. Future studies will focus on investigating N_d-LWP sensitivity and its effect on aerosol perturbation using Lagrangian cloud tracking, which can improve the understanding of the aerosol effect on LWP sensitivity.

Appendix A: Derivation of the adiabatic liquid water path, LWP_ad

Following Toledo et al. (2021), the adiabatic liquid water path is computed from the standard adiabatic assumption, in which the adiabatic LWC increases linearly with height above cloud-base and corresponding equation is given by,

$${\mathrm{LWP}}_{\mathrm{ad}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\gamma }}_{\mathrm{ad}}{H}^{\mathrm{2}},$$

where γ_ad is the vertical gradient of adiabatic LWC, and the cloud depth H is

$$H=z_t-z_{\mathrm{b}}{z}_{\mathrm{b}}\le z\le z_t\left(m\right)$$

The adiabatic LWC gradient γ_ad is evaluated at cloud-base using the cloud-base temperature T_b (K) and the pressure p_b (Pa):

$${\mathit{\gamma }}_{\mathrm{ad}}={\mathit{\rho }}_{\mathrm{d}}\left[\left({\displaystyle \frac{(\mathit{\epsilon }+w_{\mathrm{s}})w_{\mathrm{s}}{L}_{\mathrm{v}}}{R_{\mathrm{d}}{T}_{\mathrm{b}}^{\mathrm{2}}}}\right){\mathrm{\Gamma }}_w-{\displaystyle \frac{g{w}_{\mathrm{s}}{p}_{\mathrm{b}}}{(p_{\mathrm{b}}-e_{\mathrm{s}})R_{\mathrm{d}}{T}_{\mathrm{b}}}}\right],$$

where Γ_w is the moist-adiabatic temperature lapse rate (units: K m⁻¹), ρ_d is the dry-air density (units: kg m⁻³), w_s is the saturation mixing ratio (units: kg kg⁻¹), e_s=e_s(T_b) is the saturation vapour pressure (units: Pa) at T_b (cloud base temperature), and p_b (cloud base pressure).

The moist-adiabatic lapse rate is given by

$${\mathrm{\Gamma }}_w={\displaystyle \frac{g}{c_{\mathrm{p}}}}{\displaystyle \frac{\mathrm{1}+{\displaystyle \frac{L_{\mathrm{v}}{w}_{\mathrm{s}}}{R_{\mathrm{d}}{T}_{\mathrm{b}}}}}{\mathrm{1}+{\displaystyle \frac{\mathit{\epsilon }{L}_{\mathrm{v}}^{\mathrm{2}}{w}_{\mathrm{s}}}{R_{\mathrm{d}}{c}_{\mathrm{p}}{T}_{\mathrm{b}}^{\mathrm{2}}}}}},$$

the dry-air density by

$${\mathit{\rho }}_{\mathrm{d}}={\displaystyle \frac{p_{\mathrm{b}}-e_{\mathrm{s}}}{R_{\mathrm{d}}{T}_{\mathrm{b}}}},$$

and the saturation mixing ratio by

$$w_{\mathrm{s}}=\mathit{\epsilon }{\displaystyle \frac{e_{\mathrm{s}}}{p_{\mathrm{b}}-e_{\mathrm{s}}}},$$

with $\mathit{\epsilon }=R_{\mathrm{d}}/R_{\mathrm{v}}$.

The saturation vapour pressure is computed using the Bolton formulation (Bolton, 1980):

$$e_{\mathrm{s}}\left(T_{\mathrm{b}}\right)=\mathrm{611.2}\exp \left({\displaystyle \frac{\mathrm{17.67}(T_{\mathrm{b}}-\mathrm{273.15})}{T_{\mathrm{b}}-\mathrm{29.65}}}\right)\mathrm{Pa}.$$

The constants used are

$$\begin{array}{rl}{\displaystyle }{R}_{\mathrm{d}}& {\displaystyle }=\mathrm{287}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}\text{(specific gas constant for dry air)},\\ {\displaystyle }{R}_{\mathrm{v}}& {\displaystyle }=\mathrm{461.5}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}\text{(specific gas constant for water vapour)},\\ {\displaystyle }g& {\displaystyle }=\mathrm{9.81}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}\text{(gravitational acceleration)},\\ {\displaystyle }{c}_{\mathrm{p}}& {\displaystyle }=\mathrm{1005}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}\text{(specific heat at constant pressure)},\\ {\displaystyle }{L}_{\mathrm{v}}& {\displaystyle }=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}\text{(latent heat of vaporization)}.\end{array}$$

Figure A1 The N_d-bin mean (a) $\stackrel{\mathrm{‾}}{\mathrm{LWC}}$ and (b) Specific humidity ($\stackrel{\mathrm{‾}}{q_{\mathrm{s}}}$) for the 2013 and 1985 simulations. The solid lines represent the smoothed mean of the mean LWC and specific humidity at certain N_d bins (same as Fig. 1). The shaded region represents the rolling standard deviation of the respective N_d-bin mean values. The green line denotes the 2013 simulation using present-day (2013) CCN concentrations, while the 1985 simulation applies CCN concentrations representative of peak aerosol loading over Europe around 1985.

Figure A2 The N_d-bin mean (a) $\stackrel{\mathrm{‾}}{\frac{\mathrm{d}T}{\mathrm{d}t}}$, (b) $\stackrel{\mathrm{‾}}{\frac{\mathrm{d}q}{\mathrm{d}t}}$, (c) advection component of temperature, and (d) advection component of moisture for the 2013 and 1985 simulations. The solid lines represent the smoothed mean of the mean values of the above variables at certain N_d bins (same as Fig. 1). The shaded region represents the rolling standard deviation of the respective N_d-bin mean values. The green line denotes the 2013 simulation using present-day (2013) CCN concentrations, while the 1985 simulation applies CCN concentrations representative of peak aerosol loading over Europe around 1985.