**Research article**
13 Aug 2019

**Research article** | 13 Aug 2019

# An emulator approach to stratocumulus susceptibility

Franziska Glassmeier Fabian Hoffmann Jill S. Johnson Takanobu Yamaguchi Ken S. Carslaw and Graham Feingold

^{1,4},

^{1,2},

^{3},

^{1,2},

^{3},

^{1}

**Franziska Glassmeier et al.**Franziska Glassmeier Fabian Hoffmann Jill S. Johnson Takanobu Yamaguchi Ken S. Carslaw and Graham Feingold

^{1,4},

^{1,2},

^{3},

^{1,2},

^{3},

^{1}

^{1}Chemical Sciences Division, NOAA Earth System Research Laboratory, 325 Broadway, Boulder, CO 80302, USA^{2}Cooperative Institute for Research in Environmental Sciences, University of Colorado Boulder, Boulder, CO 80309, USA^{3}School of Earth and Environment, University of Leeds, Woodhouse Lane Leeds, LS2 9JT, UK^{4}Department of Environmental Sciences, Wageningen University, P.O. Box 47, 6700AA Wageningen, the Netherlands

^{1}Chemical Sciences Division, NOAA Earth System Research Laboratory, 325 Broadway, Boulder, CO 80302, USA^{2}Cooperative Institute for Research in Environmental Sciences, University of Colorado Boulder, Boulder, CO 80309, USA^{3}School of Earth and Environment, University of Leeds, Woodhouse Lane Leeds, LS2 9JT, UK^{4}Department of Environmental Sciences, Wageningen University, P.O. Box 47, 6700AA Wageningen, the Netherlands

**Correspondence**: Franziska Glassmeier (franziska.glassmeier@noaa.gov)

**Correspondence**: Franziska Glassmeier (franziska.glassmeier@noaa.gov)

Received: 24 Dec 2018 – Discussion started: 21 Jan 2019 – Revised: 17 May 2019 – Accepted: 04 Jun 2019 – Published: 13 Aug 2019

The climatic relevance of aerosol–cloud interactions depends on the sensitivity of the radiative effect of clouds to cloud droplet number *N*, and liquid water path LWP.
We derive the dependence of cloud fraction CF, cloud albedo *A*_{C}, and the relative cloud radiative effect $\text{rCRE}=\text{CF}\cdot {A}_{\text{C}}$ on *N* and LWP from 159 large-eddy simulations of nocturnal stratocumulus.
These simulations vary in their initial conditions for temperature, moisture, boundary-layer height, and aerosol concentration but share boundary conditions for surface fluxes and subsidence.
Our approach is based on Gaussian-process emulation, a statistical technique related to machine learning.
We succeed in building emulators that accurately predict simulated values of CF, *A*_{C}, and rCRE for given values of *N* and LWP.
Emulator-derived susceptibilities $\partial \mathrm{ln}\text{rCRE}/\partial \mathrm{ln}N$ and $\partial \mathrm{ln}\text{rCRE}/\partial \mathrm{ln}\text{LWP}$ cover the nondrizzling, fully overcast regime as well as the drizzling regime with broken cloud cover.
Theoretical results, which are limited to the nondrizzling regime, are reproduced.
The susceptibility $\partial \mathrm{ln}\text{rCRE}/\partial \mathrm{ln}N$ captures the strong sensitivity of the cloud radiative effect to cloud fraction, while the susceptibility $\partial \mathrm{ln}\text{rCRE}/\partial \mathrm{ln}\text{LWP}$ describes the influence of cloud amount on cloud albedo irrespective of cloud fraction.
Our emulation-based approach provides a powerful tool for summarizing complex data in a simple framework that captures the sensitivities of cloud-field properties over a wide range of states.

Aerosol perturbations can lead to changes in cloud brightness and amount via the influence of aerosol on cloud formation and various aerosol–cloud interaction (ACI) processes. Our process understanding of ACI has improved greatly over recent decades; however, the radiative forcing due to ACI continues to dominate the uncertainty margin of the total anthropogenic forcing of the climate system. Due to their abundance and location, forcing and forcing uncertainty are dominated by shallow, warm clouds in marine boundary layers (Myhre et al., 2013; Boucher et al., 2013).

ACIs are notoriously hard to quantify because they pose a multiscale problem, not only in terms of spatial and temporal scales but also in terms of the effective degrees of freedom used to describe ACI in different settings.
The multiscale spectrum of approaches to the ACI problem has been described to range from “Darwinian” (low-level, high-dimensional, complex, *reductionist*) to “Newtonian” (high-level, low-dimensional, effective, *emergent*) descriptions (Harte, 2002; Feingold et al., 2016; Mülmenstädt and Feingold, 2018).

We illustrate this notion by discussing the *relative cloud radiative effect* (rCRE) as quantified by

where *F* denotes downwelling SW radiative fluxes at the surface under clear-sky (index clr) and all-sky (index all) conditions and *A*_{C} and CF denote cloud albedo and cloud fraction, respectively (Xie and Liu, 2013).
By describing their effect on the radiation budget, rCRE captures the effect of clouds on climate.
In the spirit of Platnick and Twomey (1994), ACI can be quantified based on the *susceptibility*, or normalized sensitivity, of the rCRE to the cloud droplet number concentration *N*:

The decomposition on the right-hand side is motivated by distinguishing cloud microstructure as captured by *N* from macrostructure as represented by the liquid water path, LWP.
Note that this decomposition does not necessarily align with the contributions of *A*_{C} and CF to rCRE.
Aerosol effects on cloud microstructure are associated with cloud brightness (Twomey, 1977, 1974; Boers and Mitchell, 1994; Feingold et al., 1997; Christensen and Stephens, 2011; McGibbon and Bretherton, 2017), while aerosol effects on the macrostructure of the cloud correspond to cloud amount (Albrecht, 1989; Matheson et al., 2005; Kaufman et al., 2005; Small et al., 2009; Stevens and Feingold, 2009; Zheng et al., 2010; Christensen and Stephens, 2011; Chen et al., 2014; Feingold et al., 2015; McGibbon and Bretherton, 2017).

The reductionist approach to deriving rCRE(LWP,*N*) and the partial derivatives in Eq. (2) starts by asking how *A*_{C}(*τ*) depends on cloud optical thickness *τ*,
followed by deriving the dependence of *τ*(LWP,*N*) on LWP and *N*,
which are, in turn, functions of the meteorological and aerosol conditions.
Even more complex chains of dependencies can be derived for CF(LWP,*N*) and LWP(*N*) because these relationships are shaped by the entire cloud field.
The advantage of this approach is that each link can in principle be deduced from detailed process understanding.
The disadvantage is that a large number of variables and processes need to be quantified.

The emergence-based alternative is to subsume process complexity in low-dimensional relationships that effectively describe rCRE as a function of a small set of controlling variables.
In other words, this means searching for a simple relationship for rCRE(LWP,*N*).
Alternatively, it may mean abandoning the idea of disentangling aerosol effects on cloud microstructure (*N*) and macrostructure (LWP) altogether
– similar to the definition of effective radiative forcing in Myhre et al. (2013) –
and quantifying dln rCRE∕dln *N* directly.
While the direct nature of this approach is an obvious advantage, the challenge of the emergence-based approach lies in its data-mining exploratory nature, and lack of a priori guidance for finding and justifying emergent relationships.

In this paper, we aim to combine the reductionist and emergence-based approaches to determine the partial derivatives in Eq. (2); the LWP adjustment dln LWP∕dln *N* will be the topic of an upcoming paper.
We demonstrate how a statistical method related to machine learning can be applied to derive system-wide relationships from the detailed process representation that is ingrained in a set of model simulations.
Our contribution addresses an increasing interest in machine learning approaches within the atmospheric sciences, especially in the context of parameterizing shallow clouds (Krasnopolsky et al., 2013; Schneider et al., 2017; Brenowitz and Bretherton, 2018; Gentine et al., 2018; O'Gorman and Dwyer, 2018).
This interest in utilizing modern computational statistical methods illustrates a community need to address a certain mismatch between traditional process-based cloud research and synthesizing approaches, especially for representing clouds in climate models.

We specifically apply the tool of *Gaussian-process emulation* (Rasmussen and Williams, 2006; O'Hagan, 2006).
Emulation is an established method used to extract multidimensional relationships from sparse data.
It can be considered a form of kernel-based supervised machine learning.
In the atmospheric sciences, emulation has so far been used to investigate the relationship between model response and uncertain parameters associated with physical parameterizations and to a lesser extent boundary conditions, e.g., Lee et al. (2011, 2013), Johnson et al. (2015), and Posselt et al. (2016).
We adapt this method to quantitatively derive relationships between cloud-field properties that evolve over the course of numerical simulations, namely rCRE(LWP,*N*), *A*_{C}(LWP,*N*), and CF(LWP,*N*).

We present state-of-the-art large-eddy simulations (LESs) of stratocumulus (Sect. 2) and
demonstrate that our approach (Sect. 3) successfully translates process understanding captured by the LES into
a quantification of the rCRE, *A*_{C}, and CF and their relationships to LWP and *N* (Sect. 4).
As an application, we then derive and discuss the partial susceptibilities in Eq. (2) (Sect. 5)
before we conclude (Sect. 6).

We perform LES with the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall, 2003).
Our domain measures 48 km×48 km at a horizontal resolution of 200 m and a vertical resolution of 10 m.
The time step is 1 s.
Simulations are nocturnal and of 12 h duration.
Sea surface temperature, subsidence, and surface fluxes are fixed to the values of Ackerman et al. (2009).
The model activates aerosol particles based on the prognosed supersaturation and simulates condensation and/or evaporation and collision-coalescence using a bin-emulating 2-moment approach (Feingold et al., 1998).
Particles are removed by collision-coalescence, scavenging and wet deposition. We assume a surface source of aerosol particles of 70 cm^{−2} s^{−1} (Yamaguchi et al., 2017; Kazil et al., 2011).
Integrated properties such as cloud water path (CWP), rain water path (RWP), and their sum, the LWP, are calculated directly from cloud and rain water mass mixing ratios.
Drop number concentration is calculated from the prognosed cloud and rain number concentrations and is usually dominated strongly by the former.
For more details on the model setup see Yamaguchi et al. (2017).

Following Feingold et al. (2016), we simulate 191 stratocumulus (Sc) cases with different initial conditions. We simultaneously vary six initial conditions of the Sc field: liquid water potential temperature ($\mathrm{284}<{\mathit{\theta}}_{\mathrm{l}}/\mathrm{K}<\mathrm{294}$), the total mixing ratio ($\mathrm{6.5}<{q}_{\mathrm{t}}/\left(\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{-\mathrm{1}}\right)<\mathrm{10.5}$), and the aerosol concentration ($\mathrm{30}<{N}_{\mathrm{a}}/{\mathrm{cm}}^{-\mathrm{3}}<\mathrm{500}$) in a mixed layer, as well as the initial height ($\mathrm{500}<{H}_{\mathrm{mix}}/\mathrm{m}<\mathrm{1300}$) of that mixed layer and the jumps ($\mathrm{6}<\mathrm{\Delta}{\mathit{\theta}}_{\mathrm{l}}/\mathrm{K}<\mathrm{10}$, $-\mathrm{10}<\mathrm{\Delta}{q}_{\mathrm{t}}/\left(\text{g}\phantom{\rule{0.125em}{0ex}}{\text{kg}}^{-\mathrm{1}}\right)<-\mathrm{6}$) at the inversion above the mixed layer.

To generate our ensemble of model simulations, we sample well-spaced combinations of these initial conditions over the six-dimensional space spanned by the ranges given above using a maximin Latin-hypercube design algorithm (Morris and Mitchell, 1995). This Latin-hypercube sampling ensures good coverage across the six-dimensional space of the initial conditions and prevents any spurious aerosol-meteorology covariability due to sampling (Feingold et al., 2016).

We create a Latin-hypercube sampling of 200 points.
Not all of these correspond to conditions for which cloud formation is expected.
Based on applying saturation adjustment as an a posteriori condition, we identify 191 suitable initial conditions.
For these initial conditions, actual LESs are performed.
From these simulations, we remove all those that do not sustain cloud.
We also remove simulations that would form rain within the first hour while collision-coalescence and sedimentation are still switched off for spin-up.
These simulations are characterized by unrealistically strong rain once collision-coalescence is allowed.
We identify such early-precipitating simulations based on the criterion that the maximum domain-averaged surface precipitation rate over the time series (0–12 h) exceeds 10 mm d^{−1}.

After this filtering, 159 of the initial 191 simulations remain. We discard 2 h of spin-up and build our analysis on hours 2 to 12, at output intervals of 10 min, which means that the time series from each simulation consists of 60 domain-averaged values. The total number of data points amounts to $\mathrm{60}\cdot \mathrm{159}=\mathrm{9540}$.

Figure 1 summarizes our dataset as a function of the domain-averaged liquid water path, LWP, which we define as the sum of cloud and rain water path, $\text{LWP}=\text{CWP}+\text{RWP}$, and the vertically averaged cloud droplet number concentration, *N*, in cloud columns with optical depth *τ*_{500 nm}>1, where the index indicates the considered wavelength.
To illustrate the system evolution, we discuss the four labeled trajectories in the figure.
Figure 2 illustrates the spatial arrangement of these trajectories.
The clouds of trajectory A deepen in response to longwave radiative cooling and attendant condensation.
In contrast, the thick clouds of trajectory B feature strong entrainment that leads to cloud thinning.
Cloud deepening and thinning approximately balance each other in a region indicated by the solid blue line in Fig. 1.
Trajectory C shows a cloud system with large droplets whose adiabatic volume-mean droplet radius at cloud top quickly reaches a critical value of about 12 µm associated with the onset of precipitation (Gerber, 1996).
This rapidly reduces the cloud droplet number (Fig. 1) and leads to the breakup of the cloud field (Fig. 2).
Trajectory D initially features droplet sizes that are too small for precipitation formation. Through cloud deepening similar to trajectory A, droplets grow until their size crosses the threshold value for precipitation formation. As for trajectory C, this means that the cloud droplet number starts to decrease.

To guide further discussion, we partition the state space into three regions: the upper right quadrant (first quadrant, Q1 in the following) is characterized by the absence of drizzle and evolution from the initial state towards decreasing LWP; the upper left quadrant (Q2) features no drizzle and increasing LWP; the lower part of the state space (Q34) features drizzle and decreasing droplet number.

We distinguish the two-dimensional *state space* spanned by LWP and *N* from the *parameter space* of the system.
Parameters are set externally and do not evolve in time.
The parameters of our simulations are its boundary conditions, especially sea surface temperature, subsidence, and surface fluxes.
Initial conditions could also be considered parameters.
Their role for the system's evolution is, however, somewhat ill-defined due to the spin-up process.
In real systems, a distinction between state variables and parameters is only approximate.
It requires a timescale separation so that slowly evolving variables can be considered fixed and parameter-like in comparison to fast-evolving variables.
While previous applications of emulators, e.g., Lee et al. (2013) and Johnson et al. (2015), have explored how the behavior of the system varies across the parameter space, we explore how it varies across the resulting state space.
Our choice of LWP and *N* as state variables is motivated by Eq. (2).
It is not a priori clear that the properties of cloud fields that arise from a six-dimensional set of initial conditions can be described as two-dimensional functions.
For such a reduction in dimensionality to be successful, it is required that multiple initial conditions in the six-dimensional space map onto individual points on the two-dimensional state space.
In hydrology, this circumstance is known as *equifinality* (Beven, 2005).
Our data do not perfectly, but only approximately, collapse onto the two-dimensional state space.
This imperfection manifests as noise in our data when presented in two dimensions.

We analyze our data based on the assumption that the relationship between the state-space variables and the cloud output variables can be modeled as a Gaussian process, and employ the technique of *Gaussian-process emulation* (Rasmussen and Williams, 2006).
The historical origins of this method lie in predicting the distribution of gold in South Africa, based on a small sample of carefully located test drillings.
Gaussian-process emulation is a preferred interpolation technique for sparse data (that is ideally well spaced) of variables that vary smoothly (no discontinuities) across the dimensions of interest.
Our data sample the LWP–*N* state space partially in a sparse manner (when considering different runs) and partially in a reasonably dense manner (within the time series of individual runs).
In contrast to the six-dimensional initial conditions, which we sampled using a space-filling Latin hypercube design, the data associated with the system evolution is not Latin-hypercube sampled.
As the system evolves and moves towards the solid blue line in Fig. 1, the coverage of the space becomes less evenly distributed, which can lead to issues of instability in parameter estimation when fitting a Gaussian-process emulator.
We therefore adapt the emulation approach to our situation.
To fulfill the methodological requirements of sparsity and sampling, we generate a set of approximately Latin-hypercube sampled subsets of the data and construct a Gaussian-process emulator for each subset (see Sect. 3.1).
Together, this set of fitted emulators forms an emulator ensemble that takes into account most of the available data (see Sect. 3.2).

## 3.1 Subsampling

To build and validate the emulators in the emulator ensemble, we split the total dataset randomly into two equally large subsets. We use one of these to build the emulator (*training dataset*) and the other to independently test the emulator (*validation dataset*). The random splitting is based on equal, unconditioned probabilities for each data point to belong to either of the two datasets. It especially does not take into account to which time series a data point belongs and which data points from the same time series may already be in the same dataset.

Gaussian-process emulation is designed for reasonably sparse and well-spaced data over the dimensions of interest.
As discussed, our data do not take this form because they consist of many densely sampled time series that are themselves sparsely distributed in state space.
We therefore subsample from the total training data in the following way:
we create a virtual set of *n*_{trn} Latin-hypercube sample points in the LWP–*N* state space and replace each point in the Latin-hypercube sample by the geometrically closest point from our dataset.
Data points are added when their normalized Euclidean distance from the Latin-hypercube sampled point does not exceed 5∕*n*_{trn}. We do not use data points twice and training and validation data are treated as completely separate such that a data point cannot be selected for both training and validation. Figure 3 illustrates two such examples of subsampling. We proceed in the same way to select a subsample of *n*_{vld}=2*n*_{trn} validation data points from the validation dataset. The Latin hypercubes underlying the samplings for the training and the validation dataset augment one another so that their combination is also Latin-hypercube sampled. We achieve the Latin-hypercube samples using the `R`

package `lhs`

(R Core Team, 2018; Carnell, 2018).

## 3.2 Ensemble emulation and averaging

By varying the random seed of the initial Latin-hypercube sampling, we create an ensemble of subsamples.
To each subsample we apply Gaussian-process emulation, constructing and validating a separate emulator model for each subsample in turn.
For a general overview of the mathematical details of the emulator model, we refer the reader to Johnson et al. (2015).
It is based on a Bayesian statistical framework where we select an underlying mean and covariance structure that is then fitted given information from the training data.
We specifically assume a linear combination of LWP and *N* as the underlying mean function and use a Matérn covariance structure (Rasmussen and Williams, 2006).
We account for noise in our data (*nugget effect*).
Our data are noisy because our simulations do not perfectly, but only approximately, collapse onto the LWP–*N* state space.
This is illustrated by the fact that individual data points in region Q34 of Fig. 1 may differ in their value of RWP/LWP, i.e., their color, from closely neighboring points (cf. Fig. 5 to see the same for rCRE instead of RWP/LWP).
Emulators are fitted using the function `km()`

in the `R`

package `DiceKriging`

(Roustant et al., 2012; R Core Team, 2018).

As the data points chosen for different Latin-hypercube samplings are not necessarily different, the emulator ensemble members that we obtain in this way are not independent.
To limit the repeated use of data points, we relate the number of ensemble members *n*_{ensbl}, or random seeds, to the number of training data points for the individual emulators *n*_{trn} such that no more than 50 % of the total available training data ${n}_{\text{trn}}^{\text{tot}}$ are used:
${n}_{\text{ensbl}}={n}_{\text{trn}}^{\text{tot}}/\left(\mathrm{2}{n}_{\text{trn}}\right)$.

We characterize an individual emulator within the emulator ensemble by its root-mean-square error (RMSE) in predicting the validation data.
We characterize the emulator ensemble as a whole by a weighted mean *μ*, where the weighting *w* depends on the RMSE of the individual emulators:

For the example of the rCRE emulator shown in Fig. 5, Fig. 3 illustrates the best and worst sampling within the ensemble as quantified by the RMSE of the corresponding individual emulators.

Figure 4 shows the spread of emulator RMSE that is obtained when varying the number of training data points *n*_{trn} used to build an emulator ensemble.
The median RMSE tends to decrease with an increasing number of data points, while the number of ensemble members decreases commensurately.
As a compromise between the quality of individual emulators and ensemble statistics we choose *n*_{trn}=50.

Figures 5 and 6 demonstrate the successful application of our emulator ensemble technique to derive surfaces of rCRE, *A*_{c}, and CF (all determined using a threshold of *τ*>1) as a function of LWP and *N* from the multi-time-series data shown in Fig. 1.
The emulated surfaces successfully predict simulation outcomes (validation data) that were not used in creating the emulator (training data) as shown by scatter plots and data points in the figures.

In accordance with Eq. (1), the shape of the emulated rCRE surface follows from the surfaces of cloud albedo and cloud fraction (Fig. 6a, c).
Cloud fraction generally decreases with decreasing LWP as the Sc deck entrains and thins until the detection threshold, *τ*=1, occurs.
Its surface is dominated by a gradual shift of its isolines from Q2 to Q34 as rain formation sets in.
In the shift region (Q3, left part of Q34), cloud fraction depends strongly on *N*.
Elsewhere, isolines run mostly in the vertical direction such that CF is largely independent of *N*.
Cloud albedo is characterized by negatively sloped isolines so that cloud albedo tends to increase with both LWP and *N*.
The isolines, and thus the dependence of cloud albedo on LWP and *N*, are distorted in the drizzling region (Q34).

This distortion in Q3 is caused by the bimodality of the drop-size distributions, so that cloud albedo and fraction are influenced by the radiative effects of cloud droplets as well as rain drops.
Figure 6b and d considers the cloud water path (CWP) instead of LWP as *x* axis.
This transformation leads to a shift in the isolines:
cloud fraction isolines become approximately vertical as *τ*>1 is controlled by CWP and hardly influenced by the additional contribution of RWP to LWP.
The tilt of the cloud albedo isolines is reversed, indicating a contribution of RWP to total cloud albedo.

## 4.1 Uncertainty

Our approach features five different kinds of uncertainties, or errors. First, the Gaussian-process emulation returns a random variable, i.e., a probability distribution of possible surfaces. This uncertainty depends on the training data used to build an individual emulator. It can be quantified by a standard deviation and its values can be inferred from Fig. 8.

Second, the quality of the emulator mean function is quantified by the RMSE of an emulator in predicting the validation data.
This measure of uncertainty depends on the training as well as the validation data that is used for a specific emulator, i.e., one member of the emulator ensemble.
It is not spatially resolved but averages the error of the emulator over the whole LWP–*N* state space.
As indicated in the caption of Fig. 3, RMSE lies between 0.01 and 0.02 for the rCRE emulators.

Third, we have the error due to the noise that arises when collapsing our dataset onto the LWP–*N* state space.
This uncertainty is the most fundamental uncertainty because it cannot be reduced by increasing the number of data points available.
It is inherent to our modeling of rCRE as two-dimensional function of LWP and *N*.

Fourth, we have the uncertainty of the emulator ensemble. Following Eq. (3), we quantify this uncertainty using a weighted ensemble standard deviation,

where the weights *w*_{i} are defined in Eq. (3) and depend on the RMSE.
Hatching in Fig. 5 shows that the ensemble uncertainty (weighted standard deviation) for values of the rCRE is mostly smaller than 0.05 and in large regions of the state space smaller than 0.01.
Comparing the different emulators, we find that the cloud-fraction emulator ensemble is more uncertain than that of the cloud albedo (Fig. 6).
This reflects that cloud albedo only depends on the local properties of the cloud field that are directly represented by LWP and *N*, while the cloud fraction is a cloud-field property.
The cloud-fraction uncertainty is mitigated by considering the combined quantity of the rCRE.

Lastly, we have a sampling uncertainty illustrated in Fig. 7.
We indicate the level of sampling uncertainty by counting the number of trajectories that sample a specific part of the state space.
To this end, we partition the LWP–*N* state space into 60×50 bins and for each bin we count how many of our 159 trajectories contain a data point within.
This uncertainty arises because our data, especially when projected onto the two-dimensional LWP–*N* state space, are noisy.
Regions of the state space that are sampled by a single trajectory are thus less reliably represented than regions where the consideration of different trajectories attenuates the noise.
Note that the nugget effect assumed for building individual emulators cannot account for this sampling uncertainty: an undersampled region does not appear noisy.

Figure 8 compares the three spatially resolved types of error for the rCRE emulator. The standard deviation of the ensemble is mostly smaller than the standard deviation of the individual emulator. This reflects the larger amount of data and information considered when building the emulator ensemble. Comparison with the sampling uncertainty indicates that the ensemble standard deviation may be overconfident in poorly sampled regions because it cannot capture the true level of noise in these regions. Therefore, because the ensemble uncertainty is generally small, we will use the sampling uncertainty to guide the interpretation of the rCRE emulator in the following.

## 4.2 Comparison to bilinear regression and effective degrees of freedom

Previous studies have determined partial susceptibilities as in Eq. (2) from binned linear (e.g., Sena et al., 2016) or bilinear regression (Jiang et al., 2010; Glassmeier and Lohmann, 2018).
We therefore add a brief comparison of our method to bilinear regression.
Figure 9 demonstrates that bilinear regression does not capture the simulation data as well as the emulator surface.
While the coefficient of determination, *r*^{2}=0.95, and RMSE=0.05 are acceptable, the regression surfaces cannot account for the tilt of isolines discussed in Sect. 4.
As a consequence, the regression surface predicts a large region of non-physical, negative, rCRE.
The reason behind the poor performance of bilinear regression is that its three free parameters are insufficient to capture the complexity of the emulated surface.
The number of degrees of freedom required to adequately capture the complexity of the simulation data can be estimated from Fig. 4.
The value of *n*_{trn} at which the emulator RMSE levels off can be interpreted as the number of degrees of freedom of the fitted surface, in our case about 50.
Performing bilinear regression for specific bins increases the total number of free parameters in the regression.
In contrast to emulation, however, this approach is limited by the requirement of a sufficient number of data points per bin.
We therefore consider emulation to be a powerful and superior alternative to binned regression studies.

While partial susceptibilities (Eq. 2) are directly obtained as the coefficients of a bilinear regression, the emulated rCRE surface requires their derivation as finite differences of the array that represents the surface. Figure 10 shows the logarithmic derivatives of the surface shown in Fig. 5.

In the upper parts of the state space (Q1 and Q2), emulator-derived susceptibilities in Fig. 10a and b compare reasonably well to theoretical results for nondrizzling conditions (black line contours), which assume a unimodal droplet-size distribution and high-cloud fraction (Boers and Mitchell, 1994; Sena et al., 2016):

In accordance with the different prefactors, rCRE is thus more susceptible to LWP than to *N*.
Susceptibilities decrease with increasing LWP, reflecting the saturation behavior of *A*_{C} for high LWP.

While Feingold et al. (1997) discuss the effect of drizzle initiation on *A*_{C} and precipitation susceptibility, the authors are not aware of studies addressing susceptibilities of rCRE, *A*_{C}, CF, or related quantities under drizzling conditions.
Our emulator approach enables us to do so.
As discussed in the context of Fig. 6, the contribution of rain water to total LWP leads to a shift of rCRE isolines and makes cloud fraction at constant LWP a function of *N*.
This creates a maximum in $\partial \mathrm{ln}\mathrm{rCRE}/\partial \mathrm{ln}N$ for fixed LWP in the vicinity of the isoline distortion (Figure 10a).
It also explains why isolines of $\partial \mathrm{ln}\mathrm{rCRE}/\partial \mathrm{ln}\text{LWP}$ are tilted in the drizzling region.
This indicates that the partial susceptibility of rCRE to *N* not only captures the radiative effects of droplet size but also accounts for cloud fraction changes, while the partial susceptibility to LWP is comparably insensitive to the latter.

We abstain from interpreting the susceptibility in the left half of Q2; due to the low sample size in this region, the observed substructure is likely not physical.

We present a new method to summarize the detailed process representation ingrained in LES into a simple picture of aerosol–cloud interactions (Eq. 2).
We have constructed ensembles of Gaussian-process emulators to extract how cloud albedo, *A*_{C}; cloud fraction, CF; and the relative cloud radiative effect, rCRE (Eq. 1), depend on the domain-averaged liquid water path, LWP, and vertically averaged cloud droplet number concentration, *N*, from a set of 159 large-eddy simulations of nocturnal stratocumulus (Sc) with different initial conditions (Fig. 1).
The initial conditions were Latin-hypercube sampled from a six-dimensional space that took into account variations in moisture and temperature profiles, including their jumps at the inversion, as well as aerosol concentration and boundary-layer height.
The emulator–ensemble approach has enabled us to accurately capture *A*_{C}, CF, and rCRE as a function of LWP and *N* over a wide range of LWP and *N* (Figs. 5 and 6).
Our results are based on an idealized set of simulations that currently do not account for varying boundary conditions like subsidence and surface fluxes.
Taking such differences into account may lead to a broader and/or denser sampling of the LWP–*N* state space.
This would extend and improve the emulated surfaces.

Emulation provides a viable and more powerful alternative to multivariate linear regression for deriving cloud-state-dependent partial susceptibilities (Eq. 2).
We demonstrate this for the partial susceptibilities of rCRE to LWP and *N* (Fig. 10).
We reproduce theoretical results for full cloud cover and monomodal droplet size distributions and extend the known relationships into the drizzling regime.
As cloud fraction remains controlled by cloud water, the contribution of rain water to total LWP leads to a strong dependence of cloud fraction on *N*, for fixed LWP.
This dependence corresponds to a strong susceptibility of rCRE to *N* in the transition region from solid to broken cloud cover.

Our results confirm the expectation that rCRE is most susceptible to microphysical perturbations in the transition region between the high- and low-cloud fraction regime of Sc.
Our new approach allows us to clarify the interpretation of Eq. (2):
the direct contribution of droplet number changes to rCRE ($\partial \mathrm{ln}\text{rCRE}/\partial \mathrm{ln}N$) captures the effect of droplet number on cloud brightness as well as the effect on cloud fraction.
This is possible because *N* controls the rain water fraction RWP/LWP at constant LWP.
The adjustment contribution, $\partial \mathrm{ln}\text{rCRE}/\partial \mathrm{ln}\text{LWP}\cdot \text{d}\mathrm{ln}\text{LWP}/\text{d}\mathrm{ln}N$, captures the effect of cloud amount on rCRE, irrespective of its distribution onto fewer or more droplets or thicker or thinner clouds at low or high cloud cover.

The methodology presented provides a powerful tool for synthesizing detailed data into a simple predictive framework.
In this paper we have demonstrated its versatility for studying the sensitivities of cloud-field properties over a wide range of states.
Subsequent work will focus on employing this emulator approach to gain a deeper understanding of LWP adjustments, dln LWP∕dln *N*.
In general, computational statistical approaches like the one discussed here have broad potential.
They enable process modelers to explore their models beyond case studies, while at the same time they empower empiricists to better account for state dependence.
This combination of approaches provides a promising avenue for improving our understanding of the uncertainties associated with the representation of shallow clouds in climate models.

Input files and the model code for reproducing the simulation data of this study are available from the authors upon request.

FG carried out the analysis. All authors contributed to developing the basic ideas, discussing the results, and preparing the article.

Ken S. Carslaw is executive editor and Graham Feingold a co-editor of ACP. Other than this, the authors declare that they have no conflict of interests.

Franziska Glassmeier acknowledges support by a National Research Council Research Associateship award at the National Oceanic and Atmospheric Administration (NOAA), by The Branco Weiss Fellowship – Society in Science, administered by the ETH Zürich, and by a Veni grant of the Dutch Research Council (NWO). Fabian Hoffmann holds a visiting fellowship of the Cooperative Institute for Research in Environmental Sciences (CIRES) at the University of Colorado Boulder, and the NOAA/Earth System Research Laboratory. Jill S. Johnson and Ken S. Carslaw were supported by the Natural Environment Research Council (NERC) under grant NE/I020059/1 (ACID-PRUF) and the UK-China Research and Innovation Partnership Fund through the Met Office Climate Science for Service Partnership (CSSP) China as part of the Newton Fund. Ken S. Carslaw is currently a Royal Society Wolfson Research Merit Award holder. This research was partially supported by the Office of Biological and Environmental Research of the U.S. Department of Energy Atmospheric System Research Program Interagency Agreement DE-SC0016275. Marat Khairoutdinov graciously provided the SAM model.

This research has been supported by the National Research Council (Research Associateship Award), the The Branco Weiss Fellowship – Society in Science, Dutch Research Council NWO (Veni), the Cooperative Institute for Research in Environmental Sciences (CIRES) at the University of Colorado Boulder and the NOAA/Earth System Research Laboratory (Visiting Fellowship), the Natural Environment Research Council (grant no. NE/I020059/1 (ACID-PRUF)), the Newton Fund/Met Office (UK-China Research and Innovation Partnership Fund), and the Royal Society (Wolfson Merit Award).

This paper was edited by Johannes Quaas and reviewed by two anonymous referees.

Ackerman, A. S., vanZanten, M. C., Stevens, B., Savic-Jovcic, V., Bretherton, C. S., Chlond, A., Golaz, J.-C., Jiang, H., Khairoutdinov, M., Krueger, S. K., Lewellen , D. C., Lock, A., Moeng, C.-H., Nakamura, K., Petters, M. D., Snider, J. R., Weinbrecht, S., and Zulauf, M.: Large-Eddy Simulations of a Drizzling, Stratocumulus-Topped Marine Boundary Layer, Mon. Weather Rev., 137, 1083–1110, https://doi.org/10.1175/2008mwr2582.1, 2009. a

Albrecht, B. A.: Aerosols, Cloud Microphysics, and Fractional Cloudiness, Science, 245, 1227–1230, https://doi.org/10.1126/science.245.4923.1227, 1989. a

Beven, K.: A manifesto for the equifinality thesis, J. Hydrol., 320, 18–36, https://doi.org/10.1016/j.jhydrol.2005.07.007, 2005. a

Boers, R. and Mitchell, R. M.: Absorption feedback in stratocumulus clouds Influence on cloud top albedo, Tellus A, 46, 229–241, https://doi.org/10.3402/tellusa.v46i3.15476, 1994. a, b, c

Boucher, O., Randall, D., Artaxo, P., Bretherton, C., Feingold, G., Forster, P., Kerminen, V.-M., Kondo, Y., Liao, H., Lohmann, U., Rasch, P., Satheesh, S., Sherwood, S., Stevens, B., and Zhang, X.: Clouds and Aerosols, in: Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to IPCC AR5, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P., Cambridge, https://doi.org/10.1017/CBO9781107415324, 2013. a

Brenowitz, N. D. and Bretherton, C. S.: Prognostic Validation of a Neural Network Unified Physics Parameterization, Geophys. Res. Lett., 45, 6289–6298, 2018. a

Carnell, R.: lhs: Latin Hypercube Samples, available at: https://CRAN.R-project.org/package=lhs (last access: 2 July 2019), 2018. a

Chen, Y.-C., Christensen, M. W., Stephens, G. L., and Seinfeld, J. H.: Satellite-based estimate of global aerosol-cloud radiative forcing by marine warm clouds, Nat. Geosci., 7, 643–646, 2014. a

Christensen, M. W. and Stephens, G. L.: Microphysical and macrophysical responses of marine stratocumulus polluted by underlying ships: evidence of cloud deepening, J. Geophys. Res., 116, D03201, https://doi.org/10.1029/2010JD014638, 2011. a, b

Feingold, G., Boers, R., Stevens, B., and Cotton, W. R.: A modeling study of the effect of drizzle on cloud optical depth and susceptibility, J. Geophys. Res., 102, 13527–13534, 1997. a, b

Feingold, G., Walko, R. L., Stevens, B., and Cotton, W. R.: Simulations of marine stratocumulus using a new microphysical parameterization scheme, Atmos. Res., 47–48, 505–528, 1998. a

Feingold, G., Koren, I., Yamaguchi, T., and Kazil, J.: On the reversibility of transitions between closed and open cellular convection, Atmos. Chem. Phys., 15, 7351–7367, https://doi.org/10.5194/acp-15-7351-2015, 2015. a

Feingold, G., McComiskey, A., Yamaguchi, T., Johnson, J. S., Carslaw, K. S., and Schmidt, K. S.: New approaches to quantifying aerosol influence on the cloud radiative effect, P. Natl. Acad. Sci. USA, 113, 5812–5819, https://doi.org/10.1073/pnas.1514035112, 2016. a, b, c

Gentine, P., Pritchard, M., Rasp, S., Reinaudi, G., and Yacalls, G.: Could Machine Learning Break the Convection Parameterization Deadlock?, Geophys. Res. Lett., 45, 5742–5751, 2018. a

Gerber, H.: Microphysics of Marine Stratocumulus Clouds with Two Drizzle Modes, J. Atmos. Sci., 53, 1649–1662, 1996. a

Glassmeier, F. and Lohmann, U.: Precipitation susceptibility and aerosol buffering of warm and mixed-phase orographic clouds in idealized simulations, J. Atmos. Sci., 75, 1173–1194, https://doi.org/10.1175/JAS-D-17-0254.1, 2018. a

Harte, J.: Toward a Synthesis of the Newtonian and Darwinian Worldviews, Phys. Today, 55, 29–34, 2002. a

Jiang, H., Feingold, G., and Sorooshian, A.: Effect of Aerosol on the Susceptibility and Efficiency of Precipitation in Warm Trade Cumulus Clouds, J. Atmos. Sci., 67, 3525–3540, https://doi.org/10.1175/2010jas3484.1, 2010. a

Johnson, J. S., Cui, Z., Lee, L. A., Gosling, J. P., Blyth, A. M., and Carslaw, K. S.: Evaluating uncertainty in convective cloud microphysics using statistical emulation, J. Adv. Model. Earth Syst., 7, 162–187, https://doi.org/10.1002/2014MS000383, 2015. a, b, c

Kaufman, Y. J., Koren, I., Remer, L. A., Rosenfeld, D., and Rudich, Y.: The effect of smoke, dust, and pollution aerosol on shallow cloud development over the Atlantic Ocean, P. Natl. Acad. Sci. USA, 102, 11207–11212, https://doi.org/10.1073/pnas.0505191102, 2005. a

Kazil, J., Wang, H., Feingold, G., Clarke, A. D., Snider, J. R., and Bandy, A. R.: Modeling chemical and aerosol processes in the transition from closed to open cells during VOCALS-REx, Atmos. Chem. Phys., 11, 7491–7514, https://doi.org/10.5194/acp-11-7491-2011, 2011. a

Khairoutdinov, M. F. and Randall, D. A.: Cloud resolving modeling of the ARM Summer 1997 IOP: Model formulation, results, uncertainties, and sensitivities, J. Atmos. Sci., 60, 607–625, https://doi.org/10.1175/1520-0469(2003)060<0607:CRMOTA>2.0.CO;2, 2003. a

Krasnopolsky, V. M., Fox-Rabinovitz, M. S., and Belochitski, A. A.: Using ensemble of neural networks to learn stochastic convection parameterizations from climate and numerical weather prediction models from data simulated by a cloud resolving model, Advances in Artificial Neural Systems, 2013, 485913, 2013. a

Lee, L. A., Carslaw, K. S., Pringle, K. J., Mann, G. W., and Spracklen, D. V.: Emulation of a complex global aerosol model to quantify sensitivity to uncertain parameters, Atmos. Chem. Phys., 11, 12253–12273, https://doi.org/10.5194/acp-11-12253-2011, 2011. a

Lee, L. A., Pringle, K. J., Reddington, C. L., Mann, G. W., Stier, P., Spracklen, D. V., Pierce, J. R., and Carslaw, K. S.: The magnitude and causes of uncertainty in global model simulations of cloud condensation nuclei, Atmos. Chem. Phys., 13, 8879–8914, https://doi.org/10.5194/acp-13-8879-2013, 2013. a, b

Matheson, M. A., Coakley, J. A., and Tahnk, W. R.: Aerosol and cloud property relationships for summertime stratiform clouds in the northeastern Atlantic from Advanced Very High Resolution Radiometer observations, J. Geophys. Res., 110, D24204, https://doi.org/10.1029/2005JD006165, 2005. a

McGibbon, J. and Bretherton, C. S.: Skill of ship-following large-eddy simulations in reproducing MAGIC observations across the northeast Pacific stratocumulus to cumulus transition region, J. Adv. Model Earth Syst., 9, 810–831, https://doi.org/10.1002/2017MS000924, 2017. a, b

Morris, M. D. and Mitchell, T. J.: Exploratory designs for computer experiments, J. Stat. Plann. Inference, 43, 381–402, https://doi.org/10.1016/0378–3758(94)00035-T, 1995. a

Mülmenstädt, J. and Feingold, G.: The radiative forcing of aerosol-cloud interactions in liquid clouds: wrestling and embracing uncertainty, Curr. Clim. Change Rep., 4, 23–40, https://doi.org/10.1007/s40641-018-0089-y, 2018. a

Myhre, G., Shindell, D., Bréon, F.-M., Collins, W., Fuglestvedt, J., Huang, J., Koch, D., Lamarque, J.-F., Lee, D., Mendoza, B., Nakajima, T., Robock, A., Stephens, G., Takemura, T., and Zhang, H.: Anthropogenic and Natural Radiative Forcing, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to IPCC AR5, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P., Cambridge, https://doi.org/10.1017/CBO9781107415324, 2013. a, b

O'Gorman, P. A. and Dwyer, J. G.: Using machine learning to parameterize moist convection: potential for modeling of climate change and extreme events, J. Adv. Model Earth Syst., 10, 2548–2563, https://doi.org/10.1029/2018MS001351, 2018. a

O'Hagan, A.: Bayesian analysis of computer code outputs: A tutorial, Reliab. Eng. Syst. Safe., 38, 1290–1300, https://doi.org/10.1016/j.ress.2005.11.025, 2006. a

Platnick, S. and Twomey, S.: Determining the Susceptibility of Cloud Albedo to Changes in Droplet Concentration with the Advanced Very High Resolution Radiometer, J. Appl. Meteor., 33, 334–347, 1994. a

Posselt, D. J., Fryxell, B., Molod, A., and Williams, B.: Quantitative Sensitivity Analysis of Physical Parameterizations for Cases of Deep Convection in the NASA GEOS-5, J. Climate, 29, 455–479, https://doi.org/10.1175/jcli-d-15-0250.1, 2016. a

R Core Team: R: A Language and Environment for Statistical Computing,, R Foundation for Statistical Computing, Vienna, Austria, available at: https://www.R-project.org/ (last access: 4 July 2019), 2018. a, b

Rasmussen, C. E. and Williams, C. K. I.: Gaussian Processes for Machine Learning, MIT press, London, 2006. a, b, c

Roustant, O., Ginsbourger, D., and Deville, Y.: DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization, J. Stat. Softw., 51, 1–55, 2012. a

Schneider, T., Lan, S., Stuart, A., and Teixeira, J.: Earth System Modeling 2.0: A Blueprint for models that learn from observations and targeted high-resolution simulations, Geophys. Res. Lett., 44, 12396–12417, https://doi.org/10.1002/2017GL076101, 2017. a

Sena, E. T., McComiskey, A., and Feingold, G.: A long-term study of aerosol–cloud interactions and their radiative effect at the Southern Great Plains using ground-based measurements, Atmos. Chem. Phys., 16, 11301–11318, https://doi.org/10.5194/acp-16-11301-2016, 2016. a, b

Small, J. D., Chuang, P. Y., Feingold, G., and Jiang, H.: Can aerosol decrease cloud lifetime?, Geophys. Res. Lett., 36, L16806, https://doi.org/10.1029/2009GL038888, 2009. a

Stevens, B. and Feingold, G.: Untangling aerosol effects on clouds and precipitation in a buffered system, Nature, 461, 607–613, https://doi.org/10.1038/nature08281, 2009. a

Twomey, S.: Pollution and the planetary albedo, Atmos. Environ., 8, 1251–1256, https://doi.org/10.1016/0004-6981(74)90004-3, 1974. a

Twomey, S.: The Influence of Pollution on the Shortwave Albedo of Clouds, J. Atmos. Sci., 34, 1149–1152, 1977. a

Xie, Y. and Liu, Y.: A new approach for simultaneously retrieving cloud albedo and cloud fraction from surface-based shortwave radiation measurements, Environ. Res. Lett., 8, 044023, https://doi.org/10.1088/1748-9326/8/4/044023, 2013. a

Yamaguchi, T., Feingold, G., and Kazil, J.: Stratocumulus to Cumulus Transition by Drizzle, J. Adv. Model. Earth Syst., 9, 2333–2349, https://doi.org/10.1002/2017ms001104, 2017. a, b

Zheng, X., Albrecht, B., Minnis, P., Ayers, K., and Jonson, H. H.: Observed aerosol and liquid water path relationships in marine stratocumulus, Geophys. Res. Lett., 37, L17803, https://doi.org/10.1029/2010GL044095, 2010. a

- Abstract
- Introduction
- Simulations
- Building ensembles of emulators
- Emulators for rCRE, cloud albedo, and cloud fraction
- Partial susceptibilities of rCRE to droplet number and LWP
- Conclusions
- Code and data availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Simulations
- Building ensembles of emulators
- Emulators for rCRE, cloud albedo, and cloud fraction
- Partial susceptibilities of rCRE to droplet number and LWP
- Conclusions
- Code and data availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References