An aerosol activation metamodel of v 1 . 2 . 0 of the pyrcel cloud parcel model : development and offline assessment for use in an aerosol – climate model

We describe an emulator of a detailed cloud parcel model which has been trained to assess droplet nucleation from a complex, multimodal aerosol size distribution simulated by a global aerosol–climate model. The emulator is constructed using a sensitivity analysis approach (polynomial chaos expansion) which reproduces the behavior of the targeted parcel model across the full range of aerosol properties and meteorology simulated by the parent climate model. An iterative technique using aerosol fields sampled from a global model is used to identify the critical aerosol size distribution parameters necessary for accurately predicting activation. Across the large parameter space used to train them, the emulators estimate cloud droplet number concentration (CDNC) with a mean relative error of 9.2 % for aerosol populations without giant cloud condensation nuclei (CCN) and 6.9 % when including them. Versus a parcel model driven by those same aerosol fields, the best-performing emulator has a mean relative error of 4.6 %, which is comparable with two commonly used activation schemes also evaluated here (which have mean relative errors of 2.9 and 6.7 %, respectively). We identify the potential for regional biases in modeled CDNC, particularly in oceanic regimes, where our bestperforming emulator tends to overpredict by 7 %, whereas the reference activation schemes range in mean relative error from −3 to 7 %. The emulators which include the effects of giant CCN are more accurate in continental regimes (mean relative error of 0.3 %) but strongly overestimate CDNC in oceanic regimes by up to 22 %, particularly in the Southern Ocean. The biases in CDNC resulting from the subjective choice of activation scheme could potentially influence the magnitude of the indirect effect diagnosed from the model incorporating it.


Introduction
Aerosols play a critical role in the climate system by interacting with radiation through several different mechanisms.Depending on their composition, aerosol particles can directly scatter or absorb incoming solar radiation, leading to a direct radiative effect and rapid changes in the energy budgets of the surface and atmosphere.Additionally, aerosol particles mediate the production of clouds by providing surface area on which water vapor may condense to form droplets.Through this second pathway, changes in the aerosol population perturb the radiative properties of clouds by altering their microstructure and life cycle, thereby impacting the planetary radiative budget.Despite decades of focused research by the scientific community, the radiative forcing produced through this second pathway, known as aerosol-cloud interactions, remains one of the largest uncertainties in understanding contemporary and future climate change on both regional and global scales (Boucher et al., 2013).
To include this second pathway, contemporary Earth system models predict cloud droplet number concentration (CDNC) by evaluating the nucleation of droplets (aerosol activation) from their simulated aerosol fields.As a result, these models can resolve aerosol-climate indirect effects which arise when anthropogenic aerosol emissions influence cloud microphysical and optical properties through impacting the CDNC burden.The interactions between aerosol particles, water vapor, and cloud droplets are often described using the conceptual model of a possibly entraining, adiabatic cloud parcel (e.g., Feingold and Heymsfield, 1992;Nenes et al., 2001;Ervens et al., 2005;Topping et al., 2013).However, it is not practical to directly include these calculations in global Published by Copernicus Publications on behalf of the European Geosciences Union.models, because of their computational complexity and because parcel theory describes a process occurring on spatial and temporal scales much finer than those resolved by climate model grids.Some efforts have sought to incorporate fine-scale information about aerosol-cloud interactions by embedding higher-resolution models within each grid cell of a global model (for example, the "multi-scale modeling framework" of Wang et al., 2011), but these are also too coarse to resolve the fine-scale motions associated with parcel theory.To resolve these issues, global models employ aerosol activation parameterizations.
A comprehensive review of the development of these parameterizations is provided by Ghan et al. (2011).In broad terms, there are two families of such schemes: look-up tables based on detailed parcel model calculations (e.g., Saleeby and Cotton, 2004;Ward et al., 2010), and physically based parameterizations (e.g., Twomey, 1959;Ghan et al., 1993;Abdul-Razzak and Ghan, 2000;Cohard and Pinty, 2000;Nenes and Seinfeld, 2003;Fountoukis and Nenes, 2005;Ming et al., 2006;Shipway and Abel, 2010).Physically based approaches often rely on empirical fits or tuning derived from parcel model calculations, but are preferable for inclusion in regional or global climate models, because they are usually applicable to an arbitrary description (modal or sectional) of the aerosol size distribution.New look-up tables must be computed whenever additional aerosol parameters are introduced; this is a problem because a look-up table's size grows exponentially as the number of parameters included increases.But look-up tables offer some advantages, including computational efficiency and the ability to account for biases in estimating CDNC, such as kinetic limitations on droplet growth (Nenes et al., 2001), the co-condensation of organic vapors (Topping et al., 2013), and competition for water vapor uptake in the presence of multiple aerosol modes (Simpson et al., 2014).Rothenberg and Wang (2016) previously used statistical emulation to conduct a sensitivity analysis of a detailed parcel model, focusing on a single, lognormal aerosol mode.Other studies have also used emulation approaches to study response surfaces of both cloud parcel models (Partridge et al., 2011) and global models (Carslaw et al., 2013), but Rothenberg and Wang further developed a framework for using their emulators as an activation parameterization, which they compared to a high-dimensional extension of traditional look-up tables.Compared to two other physically based schemes, the emulator-activation parameterizations tended to predict CDNC more accurately versus a reference parcel model, particularly in regimes with weak updraft speeds and high aerosol number concentration, where the traditional parameterizations performed the most poorly.
In this present work, we extend this approach to develop a set of metamodels trained for the aerosol and meteorology parameter space simulated by a global aerosol-climate model.The resulting metamodels can thus be directly used as activation parameterizations inside a global climate to pre-dict online CDNC, given information about the aerosol size distributions and subgrid-scale meteorology in each model grid box.To refine the parameter space used for emulation, we present an analysis of how the size distribution parameters of each aerosol mode in our model contribute to activation dynamics and droplet nucleation.We then evaluate the performance of our emulator parameterizations versus two physically based schemes which are used in the vast majority of contemporary global models (see Table 3 of Ghan et al., 2011) and assess the impact of including each of these schemes on the computational expense of our global model.

Parent aerosol-climate model
We seek to derive aerosol activation emulators for the multimode, two-moment, mixing-state-resolving Model of Aerosols for Research of Climate (MARC; version 1.0.1 here) (Kim et al., 2008(Kim et al., , 2014)).MARC builds on the NCAR Community Earth System Model (CESM; version 1.2.2 here), which is a fully coupled global climate model with subcomponents for simulating climate processes in the land, ocean, atmosphere, and sea ice domains.In particular, MARC replaces the default modal aerosol treatment (Liu et al., 2012) in the atmosphere component of the CESM, the Community Atmosphere Model (CAM5; version 5.3 here; Neale et al., 2012), with a scheme which simultaneously resolves both an external mixture of different aerosol species and internal mixtures between others (Wilson et al., 2001).
Table 1 summarizes the aerosol modes predicted by MARC, which includes trimodal sulfate (nucleation -NUC, Aitken -AIT, and accumulation -ACC), discrete modes for pure black carbon (BC) and a generic organic carbon species (OC), and two "mixed" modes, one of each for mixed sulfate -black carbon (MBS) and sulfate -organic carbon particles (MOS).For each of these seven modes, MARC predicts total particle mass (M) and number concentrations (N ) for a corresponding lognormal size distribution with a prescribed width (geometric standard deviation; σ g ), as well as the total mass of carbon in each of the MOS and MBS modes.Additionally, sea salt (SSLTn) and dust (DSTn) particle size distributions are each described within MARC using a four-bin, single-moment scheme with fixed particle sizes and an assumed σ g .Each mode has a prescribed particle hygroscopicity following the κ-Köhler theory (Petters and Kreidenweis, 2007) except for the MOS mode, which has a compositiondependent κ computed assuming that the carbon and sulfate in the particle form a simple internal mixture.Unlike the MOS mode, the MBS mode assumes a core-shell structure with sulfate coating a black carbon nucleus and has a fixed hygroscopicity corresponding to that of sulfate, κ = 0.507.The sea salt modes are assumed to be comprised of NaCl with κ = 1.16, and the dust modes are assumed to be a mix-ture of minerals with a hygroscopicity of κ = 0.14 (Scanza et al., 2014).The OC and BC modes are assumed to be nonhygroscopic and not significant players in aerosol activation.
The aerosol size distributions predicted by MARC interact with both radiation and cloud microphysics.MARC uses a two-moment, stratiform cloud microphysics scheme (Morrison et al., 2008) which includes an explicit source of cloud droplet formation via aerosol activation.This is facilitated by means of a physical parameterization which takes as input both the physical and chemical properties of the ambient aerosol as well as limited information about meteorology, such as the distribution of subgrid-scale vertical velocities.In CAM5 (as used by MARC), a single characteristic updraft velocity (V ), diagnosed from the grid cell turbulent kinetic energy (TKE) provided by the moist turbulence scheme (Park and Bretherton, 2009) and assumed to be isotropic, is used to estimate droplet nucleation following Ghan et al. (1997) and Lohmann et al. (1999), such that where V is the large-scale resolved updraft velocity.Furthermore, we limit 0.2 m s −1 < V < 10 m s −1 because the processes driving turbulence are not represented well in MARC, particularly those driven by cloud-top radiative cooling above the planetary boundary layer (Ghan et al., 1997).Morales and Nenes (2010) and West et al. (2014) have explored how using different characteristic updraft velocities to represent subgrid-scale variability can influence simulated aerosol-cloud interactions; in particular, West et al. (2014) showed that using a similar TKE-based parameterization produced more realistic spatial and temporal variability in V , but tends to produce an unrealistically high frequency of its minimum-permissible value.We further assume that activation occurs in non-entraining, adiabatic updrafts which carry air up and into the base of stratiform clouds.

Aerosol activation parameter space
MARC requires 24 parameters to fully describe its 15-mode aerosol size distribution, although not all of these modes are important for aerosol activation.Activation calculations additionally require three meteorological parameters (temperature, pressure, and V ).To assess this parameter space, we sample 70 snapshots of instantaneous aerosol and meteorology fields produced by a simulation of MARC forced with present-day aerosol and precursor gas emissions (Kim et al., 2014), and with interactive dust (Albani et al., 2014) and sea salt (Mårtensson et al., 2003) emissions.Only grid cells located below 500 mb and between 70 • S and 70 • N were included in this assessment, since this is where liquid clouds and droplet nucleation are most prevalent in the model.As shown in Fig. 1, V tends to take the lower bound of 0.2 m s −1 nearly 50 % of the time in this sampling dataset and rarely exceeds 1 m s −1 .On average, in continental regimes, V is slightly larger (0.41 m s −1 versus 0.32 m s −1 in ocean regimes) but has a long, positive tail maxing out between 3 and 4 m s −1 .The upper bound on V is never reached in our model output sampling.
The distributions of key aerosol size distribution parameters are shown in Fig. 2. Because particle size is a critical parameter in assessing aerosol activation, we show distributions of µ g , the geometric mean radius, instead of M.Many size distribution parameters can vary over several orders of magnitude across the globe in even a single time step.Values for N and µ g in each mode tend to covary, though there is significant heterogeneity in the overall aerosol mixing state (not shown here).For most modes, the parameter distributions have long tails where N and µ g are both small.In these cases, few particles are likely to activate due to their small size and would yield few cloud droplets even if they did.The number concentrations of DST-and SSLT-mode particles are typically very small compared to those in the ACC, MOS, and MBS modes.Aitken-mode particles in MARC are generally small but numerous, based on how the mode is defined in the model (Table 1).We note that, on average, MOS-mode particles have an intermediate κ of 0.27 ± 0.04 (standard deviation) which varies with local sulfate abundance.
The emulation method used by Rothenberg and Wang (2016) is designed to work with an arbitrary input parameter space.However, we restrict the aerosol activation emulation parameter space by eliminating parameters which exert little or no influence on the activation process.This reduces the number of parcel model simulations necessary to train emulators as well as their complexity.The OC and BC modes are assumed to be hydrophobic and hence cannot serve as cloud condensation nuclei (CCN); NUC particles are typically too small to serve as CCN, although they could be present in a large enough number to influence supersaturation and droplet nucleation.However, nucleation-mode particles are usually most abundant when larger sulfate particles (especially ACC) are also present in large concentrations and have a greater inwww.geosci-model-dev.net/10/1817/2017/Geosci.Model Dev., 10, 1817-1833, 2017  fluence on activation dynamics.We neglect size distribution parameters from these modes to build our emulators.
To further reduce the emulation parameter space, we assess the relative importance of each individual aerosol mode and its influence on activation dynamics using an ensemble of iterative, single-mode activation calculations using a detailed reference parcel model (Rothenberg and Wang, 2016).For each of our n aerosol modes, we pick a test mode and run a parcel model simulation to compute the supersaturation maximum (S max ) achieved in an updraft in which that mode is embedded.The mode which produces the minimum S max is said to be the "dominant" mode, and we record its size distribution.We then remove that aerosol from the original set of n modes and revisit each of the n − 1 remaining modes, including now the first dominant mode along with each new test mode.The end result of n − 1 iterations is a sorting of the modes based on their contribution to reducing S max in the parcel model simulations.
We apply this algorithm to a set of 50 000 aerosol size distribution and meteorology parameters taken from our reference MARC simulation.Overall, ACC is the dominant mode in 96.5 % of the sample cases.Infrequently, MOS and small dust (DST1) are the dominant modes, accounting for all of the remaining cases.When ACC dominates the activation dynamics, MBS/MOS or the smallest sea salt mode (SSLT1) is the second most dominant (in 10.3, 36.2, and 52.8 % of cases, respectively).In 85 % of all the sample cases, three of ACC, MOS, MBS, or SSLT01 comprise the top three dominant modes.
Figure 3 illustrates the potential error in calculating both S max and N act (number concentration of activated particles or CDNC) for each iteration of the activation calculations relative to using the complete MARC aerosol population, aggregated by the first dominant mode.Using a subset of the modes leads to an overprediction of N act due to predicting a high value for S max .This is consistent with the physics of the activation problem; the presence of more aerosol surface area on which condensation can occur tends to deplete water vapor in the ascending parcel more quickly, suppressing the development of a higher S max .But as Fig. 3 shows, this overprediction decreases rapidly as additional modes are included in the calculation.Part of this decrease is related to the fact that adding modes in each iteration captures a higher fraction of the total aerosol number; on average, the first dominant mode contains 70 % ± 27 % of the total aerosol number, which increases to 80 %±17 % and 89 %±13 % after adding the second and third modes.Following this increase in fraction of the total aerosol number included by the dominant mode set in each iteration is a decrease in the absolute error in N act relative to the full aerosol population, with an average of less than 1 cm −3 and maximum of 57 cm −3 by the third iteration.Although giant CCN particles are well represented in MARC by dust and sea salt modes, they are typically too few in number to largely influence droplet nucleation, except in rare cases where other modes are not abundant.
Based on the sampling and iterative calculations presented in this section, we define the activation emulation parameter space as in Table 2.We consider the full range of potential values for V but further restrict the emulation to include the most frequent dominant aerosol modes.Table 2 further defines ranges considered for each of the aerosol and meteorological parameters, and indicates the percentile of the sampling data at which the extreme values of each range occur.In all cases, we include upper bounds above the 98th percentile of the sampling data.For the giant CCN modes (DST and SSLT), we restrict the lower ranges of N but at values where small changes in N could only have very small contributions to changes in N act .

Emulator construction
The following sections briefly describe the chosen cloud parcel model and emulation technique, polynomial chaos expansion.For more details on both techniques and their application, we refer the reader to Rothenberg and Wang (2016), which derives activation emulators for a simplified, single lognormal aerosol mode.

Parcel model
Adiabatic cloud parcel models are a standard modeling tool for detailed assessments of aerosol activation and other studies focused on the composition of atmospheric particulates (Seinfeld and Pandis, 2006).In such a model, a constantspeed updraft drives adiabatic cooling in a closed, zerodimensional air parcel within which is any number or configuration of aerosol particles.Initially prescribed a temperature, pressure, and water vapor content, the cooling parcel eventually develops a supersaturation with respect to water vapor.In a sufficiently supersaturated environment, water vapor condenses on particulate surfaces.This condensation releases latent heat and slows down the cooling of the parcel but more importantly acts to sink water vapor mass.In terms of producing supersaturation, these counteracting processes can be expressed as where V is the updraft speed, α(T , P ) = (gM w L/c p RT 2 ) − (gM a /RT ) and γ (T , P ) are functions weakly dependent on temperature and pressure (Leaitch et al., 1986), M w and M a are the molecular weights of water and air, L is the latent heat of vaporization of water, c p is the specific heat of dry air at constant pressure, R is the universal gas constant, g is the acceleration due to gravity, e s is the saturation vapor pressure, and w c is the liquid cloud water mass mixing ratio (please refer to Appendix A of Rothenberg and Wang, 2016, for more details).At some time t, latent heat release and cooling due to the parcel's adiabatic ascent will approximately balance such that dS dt = 0 and a supersaturation maximum, S max , will occur.Thereafter, S generally decreases, relaxing to some value close to unity as condensation drives droplet growth, quenching the ambient water vapor surplus.Beyond this point, the aerosol bifurcates into two populations: new cloud droplets which will continue to grow due to condensation and eventually collision and coalescence, and interstitial haze particles which may have taken up some water according to their hygroscopic mass, but upon which further condensation is not thermodynamically favorable.

Polynomial chaos expansion
We emulate the behavior of the detailed parcel model by applying the probabilistic collocation method (PCM; Tatang et al., 1997).PCM is a method of polynomial chaos expansion which seeks to construct a model response surface by mapping input parameters related to the initial conditions and behavior of a model to some response measured from the model.This process yields a computationally efficient yet accurate reproduction of the model.
The PCM is a non-intrusive technique which does not require modifications to an existing model in order to be applied.Instead, the PCM treats the original, full-complexity model as a black box and the chosen set of M input parameters as independent, random variables, X = X 1 , . .., X M , each with an associated probability density function (PDF).This PDF is used as a weighting function to derive a family of orthogonal polynomials which are used as the bases for the polynomial chaos expansion to be constructed, φ.Using a finite number of these bases and choosing some model response, R, we write the polynomial chaos expansion as (2) Such an expression has N t = P + 1 = (M + p)!/(M !p!) total terms, since a given chaos expansion of order p will contain p + 1 basis terms for each input parameter and combinations thereof.The coefficients α j are computed by evaluating the original model at a set of particular set of sample points, recording the response of the model, and solving a regression problem.Those sample points are generated by taking the roots of the orthogonal polynomials associated with each of the input parameters and their random variables.
In order to compute the polynomial chaos expansions, we use the Design Analysis Kit for Optimization and Terascale Applications (DAKOTA; Adams et al., 2014), version 6.1.This software automates the process of generating input parameter sets, sampling the full-complexity model to be emulated, and constructing the polynomial chaos expansion.Furthermore, it provides many useful statistical properties of the sample dataset and the chaos expansions themselves.

Application to MARC aerosol
We extend the idealized calculations in Rothenberg and Wang (2016) by applying the parcel model and chaos expansion technique previously described to construct emulators of aerosol activation suitable for use in MARC.Using the parameter set and value ranges summarized in Table 2, we consider two different cases for emulation: a "main" case which includes just the ACC, MBS, and MOS modes, and a giant CCN (gCCN) scheme which adds in dust and sea salt modes.We treat these two cases separately because of the nonlinear influence of gCCN on activation dynamics and cloud microphysics, which critically depend on the ambient CCN burden (Feingold et al., 1999).
For emulation, all the aerosol size distribution parameters are transformed using a logarithm, since they can take on values that span several orders of magnitude.We then construct uniform distributions with the associated ranges of values for each transformed parameter (Table 1).Using these uniform distributions, we cast all of the input parameters as independent random variables with a uniform probability density function covering the range of each corresponding distribution.
This methodology represents a compromise between using a high-fidelity representation of the aerosol size distribution parameters for our emulation, and the desire to build an emulator that can later be used in a GCM.However, we note that the distributions of aerosol parameters simulated by MARC are neither normal nor independent from one another.For instance, over remote maritime regions, total aerosol number concentration tends to be small but dominated by sea salt and small sulfate particles.In contrast, continental regions with anthropogenic emissions may feature much higher burdens of carbonaceous aerosol.Using both numerically generated orthogonal polynomials and statistical transformations, both of these complications can be handled directly (Isukapalli, 1999;Feinberg and Langtangen, 2015) but not with-out a trade-off.First, regions of the multidimensional input parameter space with low but nonzero probability of occurrence are unlikely to be sampled when building the emulator.For a sufficiently linear response function, this may not be an issue.More importantly, it becomes difficult to serialize the resulting chaos expansion in a form that can be reproduced for use later.This requires storing both the weights (coefficients) of the numerically generated orthogonal polynomials from the chaos expansion and some representation of how to reconstruct the transformed joint probability distribution formed from the input parameter probability distributions.Employing uniform probability distributions for each input parameter and assuming independence solves both of these problems, since canonical orthogonal polynomials with wellknown techniques to efficiently generate them can be used, and few additional, costly transformations of the parameter samples are necessary to run the emulator.We study the impact of this choice by evaluating samples of input parameter sets drawn from their true joint probability distribution (as simulated by MARC) in Sect.3.
These parameters are used to drive parcel model simulations where we record log 10 S max as the response variable.This value can then be used to diagnose the number concentration of droplets nucleated by assuming that any particles which experience their Köhler theory-predicted critical supersaturation activate.We note that although this does not resolve the issue of kinetic limitations on droplet growth and its potential to cause an underprediction in droplet number (Nenes et al., 2001), unlike existing activation schemes, our emulator accounts for the feedback of these effects on S max , so one avenue of its impact is lessened.
From a prediction of the S max achieved in an ascending parcel with the conditions passed to the emulator, we can then diagnose aerosol activation by rewriting the lognormal size distribution for each mode as a function of critical supersaturation (Ghan et al., 2011) to yield the expression where S m,i is the critical supersaturation for a particle of radius µ g,i from mode i.

Evaluation of emulators
We evaluate our emulators by applying them to both a synthetic sample of input parameters as well as real samples www.geosci-model-dev.net/10/1817/2017/Geosci.Model Dev., 10, 1817-1833, 2017 drawn from a MARC simulation.In all of our comparisons, we study third-and fourth-order chaos expansions both excluding (main) and including (gCCN) the coarse dust and sea salt modes.
As a reference, we compute activation statistics for each sample using both a detailed parcel model and two widely used activation schemes.The first parameterization, by Abdul-Razzak and Ghan (2000) (ARG), uses a pseudoanalytical solution to an integro-differential equation derived from the original adiabatic parcel model system; the second, by Morales Betancourt and Nenes (2014b) (MBN), applies an iterative scheme to partition the aerosol population into two subsets and uses different limits on the underlying analytical formulas to derive a maximum supersaturation.Because it requires a sequence of iterations to run, the MBN scheme is more computationally expensive than the ARG scheme but has the potential to include more detailed links between particle composition and condensation (Kumar et al., 2009) or entrainment into the parcel (Barahona and Nenes, 2007).Both the ARG and MBN schemes rely on some parameters fit to parcel model simulations conceptually similar (but different in implementation) to the one emulated here.

Input parameter space sampling
Using the parameter space defined in Table 2, n = 10 000 sample parameter sets were drawn using maximin Latin hypercube sampling (LHS).This randomized sampling method helps to ensure that the full aerosol and meteorology parameter space is studied while assessing its performance.
Figure 4 compares the performance of each emulator and the two reference activation schemes against parcel model simulations using all of the LHS samples for the main aerosol parameter sets.In the simulations, higher updraft speeds are nearly always associated with a much higher supersaturation maximum.For the emulators, accuracy tends to increase on average going from the third-order (Fig. 4a) to the fourthorder (Fig. 4b) scheme, although there is higher variance in the relative error compared to the parcel model at higher updraft speeds.With respect to the driving updraft speed, though, there is no consistent mode of bias -on average, the relative error is very low.The same does not hold true for the two reference schemes.The ARG scheme (Fig. 4c) tends to inaccurately predict supersaturation maxima at higher updraft speeds but is relatively well calibrated at lower updraft speeds, yielding a lower supersaturation maximum.On the other hand, the MBN scheme (Fig. 4d) is more accurate and better calibrated than either of the emulators or the ARG scheme, especially at higher updraft speeds, but tends to spuriously overestimate S max for weak updraft speeds.
Figure 5 extends this evaluated to diagnosed CDNC (N act ) nucleated for each predicted S max .For all the schemes, there can be substantial differences between the parcel model and each parameterization.This is particularly the case in regimes which produce smaller total CDNC, either due to a lower driving updraft velocity or a lower total aerosol number available to activate.The MBN scheme tends to consistently nucleate a higher number of droplets with respect to the parcel model, especially in situations which should have very few droplets, such as those where the total aerosol concentration is below 10 cm −3 .Although it does not have as consistent of a bias, the ARG scheme can both egregiously overpredict and underpredict CDNC, with these biases exaggerated at lower updraft speeds.By comparison, the emulators show much less overall bias.The mean error for the emulators follows that of S max and is small, but there is variance which tends to impart a small low or high bias on its estimates.
Both of these sets of plots are repeated in Figs. 6 and 7, but for the gCCN experiment.Qualitatively, the results for all parameterizations are very similar, with the same overall biases -especially for the ARG and MBN parameterizations.The emulators tend not to perform as well overall in the gCCN cases, although they are still the most highly calibrated scheme and do not have the velocity-regime errors that the MBN scheme has.In both the main and gCCN parameter sets, the MBN scheme tends to more regularly predict too many cloud droplets, save for polluted regimes giving rise to 100 cm −3 droplets where that bias reverses and the scheme has a tendency to underpredict droplet number.Neither the emulators nor the ARG show this same tendency in bias.
These differences in bias are most likely related to the choice of parcel model used in testing and building the ARG and MBN schemes; because each scheme relies on some empirical tuning to parcel model calculations, details in the implementation of each parcel model which influence its sensitivity should show up in ensemble evaluations of each activation scheme.The gCCN case is more taxing to simulate with parcel models using a Lagrangian description of the particle size distribution, because condensational growth is com-puted for each particle bin simultaneously.The stiffness ratio in this case will be extremely large, as the liquid water uptake by small particles in the main aerosol modes is much slower than those in the giant CCN modes.Although modern ODE solvers can automatically handle these scenarios, the subjective choice of which particular solver and how to discretize the giant CCN population (how many bins per mode) could influence the sensitivity of S max to changes in the model inputs and account for the differences observed here.
To better summarize the results in Figs.4-7, summary statistics on the error of each scheme versus their corresponding parcel model calculations are shown in Table 3.In both sampling cases, all of the parameterizations show a strong linear correlation (r 2 ) between their predictions and the result of the parcel model.The emulators (PCM order p) predict S max with lower overall absolute and relative error but with a much higher variance (not shown here).However, that lower error does not always translate into the emulators being the most accurate absolute predictors of N act .For the gCCN parameters, the ARG scheme predicts N act with a lower average mean relative error.In both parameter sets, the MBN scheme is the least accurate compared to the parcel model used in these sampling calculations.

MARC aerosol sampling
Although the sampling in the previous section fully exhausts the input parameter space over which aerosol activation may need to be assessed, it undoubtedly samples from aerosol and meteorological conditions which may not be likely to occur in the real world.To better understand the performance and potential bias of the emulators developed here and the www.geosci-model-dev.net/10/1817/2017/Geosci.Model Dev., 10, 1817-1833, 2017  existing activation schemes, we also studied a sample of n = 10 000 aerosol and meteorology parameter sets drawn directly from the MARC simulation previously used to study the distributions of simulated aerosol parameters, in contrast with the previous LHS sample.Error statistics from this sample, since they use model-produced parameter sets, are more representative of the real-world performance of the emulators and physically based activation schemes.All of the schemes were evaluated again using these parameter sets and the same detailed parcel model.This includes the main and gCCN emulators, which allows us to identify the importance of includ-ing the dust and sea salt modes as predictors in the chaos expansions.The parameters in these sets occasionally include values outside the ranges defined in Table 2 and studied in the previous section.These cases are more frequently associated with very low total aerosol number concentration, especially over the ocean where anthropogenic aerosols are limited and natural aerosols -which have a lower overall number burden -dominate.Because the distributions of the parameters of aerosol samples from oceanic and continental grid cells differ in this fundamental way, we break down the following analysis to reflect those differences.As in the previous sam- pling experiment, summary statistics on the performance of each emulator, alongside the ARG and MBN schemes, are detailed in Table 4. Qualitatively, all of the activation schemes perform similarly when evaluated against the MARC parameters as compared to the more generic sampling in the previous section.Figure 8 summarizes distributions of relative error in N act over land and ocean for each scheme.Neither the ARG nor the MBN scheme show much difference in error for the two regimes, although, on average, the MBN tends to underpredict N act .This underprediction usually occurs in regimes with higher updraft speeds and thus higher overall droplet number concentrations.In conditions with weaker updraft speeds, the MBN scheme instead tends to slightly overpredict N act .The ARG scheme is well calibrated in both regimes.
The emulators derived here do not fare as well as the physically based parameterizations when using the MARC samples.Both third-order schemes tend to overpredict droplet number over oceans, and underpredict it over land, but with an extremely large variance extending to ±100 %.However, including the effects of giant CCN measurably improves the performance of the third-order emulators in oceanic regimes.Increasing the order of the emulator also has a significant impact on their accuracy; the fourth-order scheme which neglects giant CCN produces smaller absolute error than the ARG and MBN schemes on average, and shows little bias between land and ocean regimes, indicating good convergence with its parent parcel model.On the other hand, the gCCN scheme has not yet converged by including fourthorder terms, even while its mean error statistics improve.Particularly troublesome is a secondary mode of extreme underprediction of droplet number of oceanic regimes, but this metric is deceptive.For very low total aerosol number concentrations, with the particle number in the single digits per cubic centimeter, the fourth-order gCCN scheme tends to predict half as many droplets as parcel model calculations indicate should form.This typically occurs when one or more of the input size distribution parameters (in particular, the number concentration) for the natural aerosol falls below the minimum threshold where the emulator was trained.When the emulators encounter inputs greater (lower) than these thresholds, they hold them to the maximum (minimum) value in their training range.This follows the assumption that the bounds for each parameter cover the entire range over which activation is sensitive to changes in that input.As a result, activation should be relatively insensitive to changes near the maximum or minimum values in the range for each parameter.With respect to number concentration, this must be the case; populations with fewer than 10 −3 cm −3 particles offer very little surface area for condensation and water vapor sink, and simply cannot exert a strong influence over the developing supersaturation in the parcel.That the fourth-order gCCN emulator produces a sensitivity that is too high in this regime suggests that statistical overfitting is occurring near the extremes of the input parameter space.
To contextualize these differences in N act bias over different geographical regimes, Fig. 9 remaps the testing samples back to the original MARC grid.Here, the difference in regional biases becomes much clearer.Virtually everywhere, the MBN scheme is biased a little low, but there is no systematic difference in this bias between land or ocean, or by geographical areas.The ARG scheme and the fourth-order main scheme show a different pattern; the ocean-land contrast in bias is clearly visible in the Northern Hemisphere.Furthermore, the bias is typically positive over maritime regions but negative over regions with anthropogenic aerosol influence.In particular, these regions include Europe and southeastern Asia -where aerosol distributions are dominated by anthropogenic sulfate and black carbon -and over north central Africa -where the aerosol is a mixture of both dust and organic carbon emissions from biomass burning.In the zonal average, the main_4 scheme is virtually identical to the ARG scheme.However, for both cases, as well as with MBN, there are larger biases over the southern parts of the oceans, where the aerosol is predominantly comprised of sea salt and smaller sulfate particles are produced indirectly through the emission of dimethyl sulfide (DMS).
Figure 9 also illustrates the poor performance of the gCCN_4 scheme, which underpredicts N act nearly everywhere, but especially so in the southern portions of the ocean basins.The consistent underprediction in this region explains the bimodal distribution over the ocean hinted at in Fig. 8.The gCCN_4 scheme does not perform too dissimilarly than the other schemes over regions with anthropogenic pollution or with mostly dust aerosol.

Implementation in MARC and computational expense
The ARG scheme is the original activation parameterization used within CAM5.3 to assess cloud droplet nucleation.We implemented the MBN scheme as an alternative in MARC, as well as an interface for chaos expansion-based schemes.
To use the emulators derived in this work, one must provide a NetCDF file which contains at least three pieces of information: an m-length vector of coefficients corresponding to each of the m terms in the expansion, a 2 × p-shape matrix which defines the lower and upper bounds for each of the input parameters, and an m×p matrix of integers correspond to the p variable expansion orders for each of the m terms in the expansion.
MARC caches these vectors and matrices in memory at startup, just as it caches several time-invariant terms used in both the ARG and the MBN schemes for each of the CCNproviding aerosol modes.
To estimate the impact of each scheme on MARC's performance, we performed a set of 3-month simulations initialized with fully spun-up aerosol and meteorology fields from a previous experiment.The simulations were conducted using 480 MPI tasks with two threads allocated to each task.Using the default configuration of MARC with the ARG scheme, the atmosphere component of the model averaged 6.1 s per model day.The MBN scheme averaged 7 % longer per model day, while the emulators tended to be comparable to the ARG scheme.Per model day, both the main schemes were comparable to within 0.4 % of the ARG scheme's performance, with the higher-order scheme costing an additional 0.16 %.Similarly, the gCCN schemes also compared similarly with the ARG scheme; the third-order scheme was 0.15 % faster than the ARG scheme, but the fourth-order scheme was 3 % slower.
Adding additional parameters to the chaos expansion underpinning the emulators would continue to add overhead to each evaluation by increasing the number of terms in the expansion.However, a larger penalty is incurred by increasing the expansion order for a given set of parameters, because this produces a much larger increase in the number of terms added to the expansion than adding a single parameter for the same order expansion.An assessment of the offline implementations of each scheme used in the analysis in the previous section yielded similar results.

Discussion and conclusions
In this work, we extend the metamodeling technique of Rothenberg and Wang (2016)  a complex, multimodal aerosol mixture simulated by a modern aerosol-climate model.Simultaneously, we characterize the performance of both our new emulators for aerosol activation and two widely used schemes from the literature, focusing on that same high-dimensional, complex aerosol parameter space.To identify the most important factors impacting activation in that complex parameter space, we apply a physically based approach to assess the sensitivity of activation statistics to the composition of the aerosol size distribution.Finally, we explore contrasts between aerosol and meteorology regimes over land and ocean, noting the potential for different biases in assessed cloud droplet number depending on the choice of activation scheme used in a particular global modeling application.
In ensembles of iterative calculations using a large sample of aerosol size distributions from a coupled aerosol-climate model, we note that, typically, a single mode tends to dominate activation or otherwise strongly predict the total number of droplets nucleated.This approach to understanding the sensitivities of activation dynamics on the underlying aerosol population is distinct from previously published approaches in the literature.For instance, Karydis et al. (2012) and Morales Betancourt and Nenes (2014a) apply an adjoint approach to derive the sensitivity of aerosol activation to perturbations in input parameters supplied to activation schemes.Detailed calculations using this approach yield a map of local sensitivities or gradients in the relationship between, for example, N act and one input parameter, while holding all others constant, and are thus difficult to interpret.The iterative calculations performed here aim instead to address the global sensitivity of activation to configurations of an aerosol population.
In terms of predicting CDNC, the accumulation-mode sulfate (ACC) alone serves as a good proxy for the activity of a full aerosol population in many cases, including in the presence of giant CCN and a wide swath of myriad updraft regimes.However, it is known that giant CCN exert a larger influence on precipitation formation in cleaner regimes (Feingold et al., 1999;Yin et al., 2000), and thus where ACC is abundant (especially over continents) the impacts of giant CCN on N act can be muted.This result is also model dependent in some sense; the ACC mode in MARC is not only ubiquitous but may inadvertently (and subjectively) take a range of mean particle sizes for which aerosol activation is especially sensitive.At the same time, the coarse dust and sea salt modes in MARC, on average, hold too small a number concentration to dramatically impact activation calculations, save for remote maritime regions far removed from anthropogenic sources.However, the presence of sea salt as one of the modes most frequently ranked in the top three influencers of activation agrees with previous results which indicate the presence of giant CCN can influence activation dynamics (e.g., Barahona et al., 2010).
The fact that a single mode can place such a strong constraint on aerosol activation is useful for attempts seeking to extend look-up table methods for building parameterizations.If two modes -an accumulation size and a coarse size -accurately predict aerosol activation, then one can con-strain the look-up table to just a few key aerosol size distribution parameters.The inclusion of variable aerosol composition would still likely make employing a look-up table in a global model unwieldy, though, necessitating more sophisticated approaches such as the metamodeling technique adopted here.
When sampling against the full training parameter space, our emulators perform capably.Neglecting the influence of the giant CCN modes, the mean relative error in predicting log 10 S max for the emulators is less than 1 %, which translates to mean relative errors for N act of 9.2 and 8.9 %.Including the giant CCN mode appears, at first, to dramatically increase the performance of the emulator, bringing those same metrics down to 0.3 and 6.9 % for the fourth-order scheme.Relative to the ARG and MBN schemes, the emulators are much more accurate on average when compared to our reference parcel model.Instead, we emphasize that the comparison of our emulators with the ARG and MBN scheme is motivated as part of a broader attempt to understand how the fundamental activation process dictates the simulated relationship between aerosol burden and CDNC in a global climate model.This relationship is critical for understanding uncertainty in the indirect effect produced by any given model, since the simulated background CDNC strongly correlates with its strength (Hoose et al., 2009).
Assessing the relative performance of activation schemes which, for all intents and purposes, perform extremely well at reproducing their own reference parcel models, is a critical step in establishing the parametric uncertainty in translating aerosol to droplet numbers and which underlies uncertainty in global model estimates of the indirect effect.
For this reason, we supplemented the evaluation of our emulators by using a second set of input parameter samples drawn from aerosol fields simulated by an aerosol-climate model.In contrast with previous studies, we use instantaneous fields in lieu of monthly or annual averages for our samples.Activation is inherently a fast process; because the microphysics schemes in aerosol-cloud models directly account for a tendency of new droplets formed via nucleation, the activation parameterization in any model will be called every time step and in every grid cell where clouds are occurring.Assessing activation schemes using temporally averaged aerosol fields risks missing some combinations of input parameters and limiting the range of values for which the scheme will need to accurately perform.
Most of the emulators and schemes tested here perform differently in oceanic and continental regimes, owing to the relative abundance of natural and anthropogenic aerosols in each.When focusing on the narrower range of aerosol parameters present in MARC (in comparison with the larger parameter space on which the emulators were trained), the emulators which explicitly account for giant CCN perform poorly, especially in maritime regimes dominated by sea salt.However, their counterpart performs nearly identical to the ARG scheme, showing a slight overprediction of N act mar-itime regimes and a slight underprediction over continents.In the global average, the emulator agrees better with the detailed parcel model than the ARG scheme.By comparison, the MBN scheme, while prone to underpredicting N act in both regimes in these calculations, shows far less variance in its prediction.This would suggest the MBN scheme actually performs extremely well -it is simply calibrated against a different baseline.Both the ARG and MBN schemes were developed using parcel models conceptually similar to the one employed here but differing in the details of their implementation.Furthermore, both schemes use some parameters computed from empirical fits to their associated parcel models.As a result, we do not expect perfect agreement between all the schemes evaluated here and instead note the importance of the relative differences between the CDNC simulated for each one.In particular, the MBN scheme does not show a large difference in relative error versus our parcel model between ocean and land regions, suggesting it is appropriately sensitive to a large range of different aerosol populations.
The results presented here have important implications for global modeling studies seeking to quantify uncertainty in the aerosol indirect effect on climate.While different activation schemes generally perform equally well when faced with idealized sets of input parameters (Ghan et al., 2011), their application in coupled aerosol-climate models may not be straightforward.Relative to parcel model calculations, activation schemes can likely show biases in predicting cloud droplet number in different regions of the world owing to spatial heterogeneity in the underlying aerosol and meteorology parameter distributions.This, in turn, could lead to biases in cloud radiative forcing and diagnosis of the indirect effect.For instance, Ghan et al. (2011) performed a pair of GCM experiments using two different schemes and observed a 10 % difference in the global average CDNC, which produced a 0.2 W m −2 difference in the indirect effect.Gantt et al. (2014) similarly showed large regional differences in CDNC when using a different set of activation schemes, leading to a spread 0.9 W m −2 in global average shortwave cloud forcing.These impacts on CDNC and indirect effect resulting from using different activation schemes will necessarily be model dependent, since the formulation of the basic activation diagnostic in each model is intertwined with regional and global variability in their simulated aerosol size distributions.
Future work should seek to systematically assess the differences in cloud microphysical processes and aerosol-cloud interactions arising from the choice of activation schemes in aerosol-climate models.As this work illustrates, employing emulators of detailed parcel model calculations which include complex chemical and physical effects on activation will aid with this task, since additional effects (e.g., changes in droplet surface tension due to organic surfactants; Abdul-Razzak and Ghan, 2004) can more easily be incorporated into parameterizations built using this method.How-

Figure 1 .
Figure 1.Distributions of model-predicted instantaneous subgridscale vertical velocity for near-surface (below 700 mb) grid cells broken down by land (red) and ocean (black) regimes.

Figure 2 .
Figure 2. Distributions of aerosol size distribution parameters for key modes simulated by MARC.Vertical dashed lines indicate the 5th and 95th percentiles of the sampling distribution for each parameter.Note that each mode name corresponds to those in Tables 1-2.

Figure 3 .
Figure3.Relative errors in S max (left) and N act (right) in subsequent iterations of the iterative activation calculations.The coloring of each box indicates which mode was the first or dominant one.In each box plot, the box encompasses the interquartile range and the whiskers extend to the 1st and 99th percentiles in the corresponding subsample.Outliers beyond these percentiles are not plotted.Labels in the legend correspond to accumulation mode (ACC), mixed organic -sulfate mode (MOS), mixed black carbon -sulfate mode (MBS), and smallest dust bin (DST01).

Figure 4 .
Figure 4. Predicted supersaturation maxima (%) from parcel model and activation parameterizations: third-order emulator (a), fourth-order emulator (b), ARG (c), and MBN (d).The main aerosol parameter set (excluding the dust and sea salt as predictor modes) was utilized here.Glyphs are shaded to denote updraft velocity corresponding to each sample draw (in m s −1 ) and are consistent for each panel.Solid black lines denote a factor-of-2 difference between predicted values from the parcel model and corresponding parameterization evaluations.

Figure 5 .
Figure 5.The same as Fig. 4 but shows the plotting of the predicted droplet number concentration (cm −3 ) nucleated for the aerosol main parameter set.

Figure 6 .
Figure 6.The same as Fig. 4 but for the gCCN parameter set.

Figure 7 .
Figure 7.The same as Fig. 5 but for the gCCN parameter set.

Figure 8 .
Figure8.Distributions of relative error in scheme prediction of N act versus detailed parcel, evaluated using samples taken from instantaneous MARC aerosol size distribution and meteorology and colored by geographical regime.The long tail of each distribution is clipped at the extrema for each scheme.The box plot in the center of each distribution shows the median and interquartile range of the total distribution of both land and ocean samples for each scheme.

Figure 9 .
Figure 9. Mean relative error in scheme prediction of N act versus detailed parcel model plotted against location on the globe where those samples originated.At each grid location, all samples across time steps and vertical levels (below 700 mb) averaged together to compute the mean.

Table 1 .
MARC aerosol-mode size distribution and composition parameters.The MOS mode (*) has a composition-dependent density and hygroscopicity which is computed using the internal mixing state of organic carbon and sulfate present at a given grid cell and time step.

Table 2 .
Input parameter space and bounds on associated uniform probability density functions used to derive polynomial chaos expansions for MARC activation.For the lower and upper bounds on the aerosol size distribution parameters, the parenthetical values denote the percentile of the distribution for that parameter at which the bound occurs.All terms are present for the main expansion; terms affixed with an ( * ) are added for the giant CCN (gCCN) expansion.

Table 3 .
Summary statistics for error in supersaturation maxima and droplet number nucleated, predicted by emulators and activation parameterization relative to corresponding simulations with a detailed parcel model.From left to right, each column represents the coefficient of determination (r 2 ), mean absolute error (MAE), mean relative error (MRE), and the normalized root mean square error (NRMSE).

Table 4 .
The same as Fig.3but for the sampling study using MARC aerosol and meteorology parameter sets.