Easy Volcanic Aerosol ( EVA v 1 . 0 ) : an idealized forcing generator for climate simulations

Stratospheric sulfate aerosols from volcanic eruptions have a significant impact on the Earth’s climate. To include the effects of volcanic eruptions in climate model simulations, the Easy Volcanic Aerosol (EVA) forcing generator provides stratospheric aerosol optical properties as a function of time, latitude, height, and wavelength for a given input list of volcanic eruption attributes. EVA is based on a parameterized three-box model of stratospheric transport and simple scaling relationships used to derive mid-visible (550 nm) aerosol optical depth and aerosol effective radius from stratospheric sulfate mass. Precalculated look-up tables computed from Mie theory are used to produce wavelength-dependent aerosol extinction, single scattering albedo, and scattering asymmetry factor values. The structural form of EVA and the tuning of its parameters are chosen to produce best agreement with the satellite-based reconstruction of stratospheric aerosol properties following the 1991 Pinatubo eruption, and with prior millennial-timescale forcing reconstructions, including the 1815 eruption of Tambora. EVA can be used to produce volcanic forcing for climate models which is based on recent observations and physical understanding but internally self-consistent over any timescale of choice. In addition, EVA is constructed so as to allow for easy modification of different aspects of aerosol properties, in order to be used in model experiments to help advance understanding of what aspects of the volcanic aerosol are important for the climate system.


Introduction
Radiative forcing by variations of stratospheric sulfate aerosol from volcanic eruptions is one of the strongest drivers of natural climate variability (Crowley, 2000;Schurer et al., 2013).To reproduce the radiative forcing of past volcanic eruptions, and thereby the related climate variability, climate model simulations require estimates of the optical properties of volcanic stratospheric aerosols.Prognostic stratospheric aerosol schemes are available; however, such schemes are computationally expensive, and many of the processes underlying them are still not well understood.For these reasons, transient simulations, such as historical or millennium simulations, usually rely on prescriptive volcanic forcing reconstructions (where "forcing" hereafter refers not specifically to "radiative forcing" but rather to any external driver of climate variability prescribed in climate model simulations).Our knowledge of past eruptions and their climate forcing is based on satellite and ground-based measurements during the recent decades, and longer-term histories can be inferred from proxies like ice cores.Different volcanic aerosol forcing sets are currently available, which use different data sources, different methodologies for combining data sources, and provide different -often incomplete -representations of aerosol properties needed for the radiative calculations of climate models.
The response of the Earth system to volcanic forcing simulated by climate models has been seen to be unrealistic in a number of prior studies.Stratospheric heating, due to the absorption of infrared radiation by volcanic aerosols, appears to be overestimated in some models (Driscoll et al., 2012).Tropospheric cooling, while relatively realistically simulated for recent eruptions (Santer et al., 2014), appears to be too Published by Copernicus Publications on behalf of the European Geosciences Union.
strong in model simulations for a number of large past eruptions, most noticeably for the eruptions of Tambora in 1815 (Brohan et al., 2012) and Samalas in 1257 (Stoffel et al., 2015).Post-volcanic anomalies of atmospheric circulation, inferred from observations, are not robustly simulated by models with prescribed forcing, either at the surface (Driscoll et al., 2012), or in the stratosphere (Charlton-Perez et al., 2013).There is also a large degree of intermodel spread in the temperature response to volcanic eruptions, and even in the radiative anomalies created by prescribed volcanic forcing (Zanchettin et al., 2016).Because there are differences in the forcing data sets used, some of which may result from differences in their implementation, it remains unclear to what degree intermodel spread in response to volcanic forcing is attributable to differences in the climate models (i.e., model uncertainty) or differences in forcing (or its implementation).
In general, to isolate model uncertainty in the response to external forcings, it is desirable to have a single forcing implementation strategy, which can be applied consistently in different models.To test the sensitivity to different aspects of the forcing, and gain understanding as to what aspects of the forcing are important for the climate response, it is further desirable to have a forcing strategy which is flexible enough to be used in sensitivity studies.Such motivations inspired the Easy Aerosol approach to prescribed tropospheric aerosols (Voigt et al., 2014), wherein the spatial structure of tropospheric aerosols was defined by simple analytical functions in latitude and longitude.Building upon the Easy Aerosol approach, the MACv2-SP module provides relatively realistic representations of anthropogenic aerosol plumes with a large degree of flexibility and utility for idealized studies (Stevens et al., 2016).
We present here a description of the Easy Volcanic Aerosol (EVA) forcing generator for use in climate model simulations.EVA provides models with the full optical properties of volcanic aerosols in terms of wavelength-dependent aerosol extinction, single scattering albedo, and asymmetry factor, given an input list of eruption locations, dates, and estimated stratospheric sulfur injections.The spatiotemporal structure of the prescribed forcing aims to strike a balance between being realistic (compared to modern observations) and generic, therefore producing consistent representations of eruptions over arbitrary time periods.The underlying parameterization is also readily modifiable, lending itself naturally to idealized sensitivity studies.EVA is comprised of a FORTRAN module that can be called directly by climate models, or can be used offline to produce forcing files which a model reads upon integration.
EVA builds directly upon the methods and results of previous volcanic aerosol forcing reconstructions, which are briefly described in Sect. 2. The EVA approach is detailed in Sect.3, and a brief comparison with other reconstructions is included in Sect. 4. A summary of EVA and outlook to potential uses and future versions is included in Sect. 5. A list of acronyms is provided in Appendix A.
2 Volcanic aerosol forcing: theory and practice Volcanic eruptions impact climate primarily though the release of sulfur gases, mostly in the form of sulfur dioxide (SO 2 ).In the atmosphere, volcanic SO 2 is chemically converted to sulfuric acid (H 2 SO 4 ), which forms liquid sulfate aerosol particles (Kremser et al., 2016).Sulfate aerosols in the stratosphere tend to have typical radii of tenths of microns (i.e., 0.1-1.0µm) (Junge et al., 1961).The size distribution of stratospheric sulfate aerosol is often approximated by the log-normal distribution, described by a distribution mean and standard deviation.From a radiative standpoint, the area-weighted mean radius, or effective radius (r eff ), is an important property of the size distribution (Hansen and Travis, 1974).
Sulfate aerosols affect the transfer of radiation through the atmosphere.Like gases or other particulate matter, sulfate aerosols can absorb and scatter incoming radiation with a spectrally dependent signature that depends on the underlying size distribution of the aerosol.Assuming spherical aerosol particles, and based on inputs describing the size distribution and the real and imaginary indices of refraction of the material, Mie theory provides an exact solution of the radiative effects of aerosols, including the proportion of incident radiation absorbed, the proportion scattered, and the variation of scattered light with direction (Hansen and Travis, 1974).There are multiple ways of representing the resulting radiative effects; a common method (e.g., Stenchikov et al., 1998) uses the following 3 parameters: aerosol extinction (EXT) represents the total attenuation of incident radiation; the single scattering albedo (SSA) represents the proportion of EXT which is scattered (as opposed to absorbed); and the scattering asymmetry factor (ASY) gives the average cosine of the scattering angle, weighted by the intensity of the scattered light as a function of the angle.It has a value of 1 for perfect forward scattering, 0 for isotropic scattering, and −1 for perfect backscatter.Finally, a very common parameter used to describe aerosol forcing is the aerosol optical depth (AOD), also called the aerosol optical thickness (AOT), which is the integral of the vertical profile of EXT.The unitless AOD therefore describes the amount by which incoming radiation is attenuated through the whole atmosphere.
The gravitational settling velocity of particles of the size of stratospheric aerosol particles is small, therefore the lifetime of sulfate aerosols in the stratosphere is significant (Junge et al., 1961).Here, "lifetime" refers to the e-folding lifetimethe characteristic timescale of a chemical or physical process in which the time derivative of a quantity is proportional to its amount -defined by the time taken for a decreasing quantity to reach 1/e of its initial amount (Jacob, 1999).Given their relatively long lifetime, the transport of sulfate aerosols is largely controlled by stratospheric dynamics.The strongest winds in the stratosphere are in the zonal (east-west) direction, which act to homogenize stratospheric composition on a fairly fast timescale (Shepherd, 2003).The Brewer-Dobson circulation (BDC) controls the distribution of stratospheric composition in the vertical-meridional plane (Holton et al., 1995;Shepherd, 2003), and is characterized by seasonally dependent two-way mixing in the midlatitudes, and a slow meridional cell of mass transport characterized by upward motion in the tropics, poleward motion in the midlatitudes, and downwelling over the poles.
The radiative effects of volcanic aerosols have been included in climate model simulations using a wide range of methods.The simplest method used involves decreasing the solar constant to reproduce the net global change in surface radiation (e.g., Metzner et al., 2014;Yoshimori et al., 2005).This method fails to reproduce the heating of the stratosphere due to aerosol absorption of infrared radiation, and also neglects spatiotemporal variation in radiative anomalies.On the other end of the spectrum, interactive stratospheric aerosol models explicitly simulate the evolution of stratospheric aerosols and their radiative impacts.While there is much to be learned from these models, significant uncertainties remain in the representation and coupling of aerosols, as evidenced by the large intermodel spread in standardized simulations of the 1815 Tambora eruption (Zanchettin et al., 2016).Furthermore, such models are computationally expensive, which usually prohibits their use for long-term (i.e., > 50-year) simulations.
The most common method of including volcanic effects in climate simulations makes use of prescribed volcanic forcing data sets.Prior works have used different methods to reconstruct volcanic forcing of climate, and implement these effects in climate model simulations.We briefly introduce reconstructions below which have been most often used in recent climate model simulations, and which are integral in the construction of EVA.

Sato/GISS
The Sato/GISS forcing data set (Sato et al., 1993(Sato et al., , 2012) ) provides stratospheric aerosol optical depth at 550 nm (AOD 550 ) and aerosol effective radius (r eff ) as a function of latitude, height, and time for the period 1850-present.It is based on a mixture of satellite observations, ground-based optical measurements, and volcanological evidence.From 2001-present, the reconstruction is based exclusively on Optical Spectrograph and InfraRed Imager System (OSIRIS) satellite measurements (Bourassa et al., 2008).For eruptions before the satellite era, the spatial structure of the AOD 550 is approximated, based on roughly scaled versions of observed eruptions, or global or hemispheric means.Aerosol effective radius is given as an empirical function of AOD 550 , based on retrievals of AOD and r eff following the Pinatubo eruption.

Ammann et al. (2003)
The Ammann et al. (2003) reconstruction provides AOD 550 as a function of latitude, height, and time for the period 1890-1999.It is based on estimates of the mass of SO 2 injected into the stratosphere, M SO 2 , from past eruptions.The spatial structure of the AOD 550 is produced by a parameterized stratospheric transport routine, representing the growth of AOD, its decay due to cross-tropopause transport, and transport from tropics to high latitudes.Compared to the Sato/GISS reconstruction, this method provides a consistent representation of volcanic forcing for all eruptions; however, it clearly simplifies the volcanic cloud evolutions compared to the observed evolutions of eruptions like Pinatubo.

Gao et al. (2008)
Building on prior work (Robock, 1981;Robock and Free, 1995), the Gao et al. (2008) reconstruction provides stratospheric sulfate aerosol mass as a function of latitude, height, and time for the period 500-2000, as well as estimates of SO 2 injection by individual volcanic events over the same period.SO 2 injections are based on ice core records of sulfate flux from Greenland and Antarctica.Using different scaling factors for tropical and extratropical eruptions based on nuclear bomb tests and modeling studies, ice-core-derived sulfate surface deposition rates are scaled to SO 2 injections (Gao et al., 2007).Sulfate aerosol mass time series, as a function of latitude and height, were produced using a modified version of the parameterized stratospheric aerosol transport scheme of Grieser and Schönwiese (1999).Conversion of sulfate aerosol mass to radiative properties was left to the implementation of the model, but a linear scaling of aerosol mass to AOD is a commonly used assumption (Schmidt et al., 2011).

Stenchikov
The Stenchikov reconstruction provides monthly mean zonal averages of stratospheric aerosol extinction, single scattering albedo, and asymmetry factor as a function of time, pressure, and wavelength, for the period from 1850 to 1999 (Driscoll et al., 2012;Schmidt et al., 2013).This data set was used by the Max Planck Institute Earth System Model and the Geophysical Fluid Dynamics Laboratory's coupled model for historical simulations as part of the Coupled Model Intercomparison Project's fifth phase (CMIP5).It is an extended version of the Pinatubo aerosol data set developed by Stenchikov et al. (1998) on the basis of satellite measurements of aerosol extinction and r eff after the Pinatubo eruption.For earlier eruptions, it is based on the Sato/GISS AOD 550 and r eff reconstruction.Stratospheric background aerosols are ignored, and only sulfate aerosols arising from volcanic eruptions are accounted for.

Crowley and Unterman (2013)
The Crowley and Unterman (2013; hereafter CU13) reconstruction, as an update to prior work (Crowley, 2000), provides AOD 550 and r eff in four equal-area latitude bands for the time period 800-2000.It is based on ice core records of sulfate flux from Greenland to Antarctica.Ice core sulfate fluxes are scaled to peak AOD 550 : for eruptions of Pinatubo magnitude and smaller, this scaling is linear and based on the ratio of satellite-observed SH AOD and the flux of sulfate to Antarctica.For stronger eruptions (like Tambora), CU13 uses nonlinear scaling first introduced by Crowley (2000), setting AOD 550 proportional to the two-thirds power of sulfate flux.AOD time series were constructed based on assumptions of linear increase to peak value, a plateau, followed by exponential decay with a 12-month timescale.As in the Sato/GISS reconstruction, effective radius is prescribed as a simple empirical function of AOD.

CCMI/SAGE_4λ
The aerosol forcing constructed for use in the Chemistry-Climate Model Initiative (CCMI; Eyring and Lamarque, 2013) provides aerosol optical properties EXT, SSA, and ASY as a function of wavelength, latitude, height, and time for the period 1960-present.It also provides internally consistent estimates of aerosol surface area density necessary for stratospheric chemistry simulations.The cornerstone of the CCMI forcing set is the four-wavelength SAGE II extinction data (SAGE_4λ), retrieval version 7, which span the period 1985-2005.During the SAGE II period, the four-wavelength extinction data are used to derive estimates of aerosol effective radius, from which the properties EXT(λ), SSA(λ) and ASY(λ) are derived assuming a single log-normal particle size distribution (Arfeuille et al., 2013).For other time periods, aerosol radiative properties are estimated mainly through single wavelength extinction retrievals from satellite instruments, and the observed correlation between mean aerosol radius and aerosol extinction during the SAGE II period.From 1979 to 1985, the reconstruction is based on single wavelength extinctions measured by the SAM II and SAGE I satellite instruments (Thomason and Peter, 2006).The 1960-1979 pre-satellite period has been constructed from SAGE-II background measurements in the late 1990s, superimposing the volcanic eruptions of Agung (1963) andFuego (1974).These eruptions were calculated by means of the AER 2-D aerosol model (Arfeuille et al., 2014;Weisenstein et al., 1997), and the results were scaled by means of stellar and solar extinction data (Stothers, 2001).The 2006-2011 period is derived from CALIPSO 532 nm backscatter data.An update to the CCMI reconstruction will be used in the historical simulations (1850-present) within the Coupled Model Intercomparison Project phase 6.
The observational data sets used in the CCMI reconstruction have data gaps, in particular when the atmosphere be-came opaque directly after volcanic eruptions, which occurred mainly in lower tropical altitudes (below 16 km), and in the high-latitude polar regions due to limited geographical sampling of the satellite instruments (Eyring and Lamarque, 2013).After the eruptions of El Chichón and Pinatubo, data gaps were filled by means of lidar ground station data and interpolation.The remaining data gaps were filled using a linear interpolation approach in altitude and latitude.
3 The EVA approach

Basis and input data
The construction of EVA relies extensively on observational constraints.Additionally, when observations cannot adequately constrain EVA parameterizations, EVA is constructed so as to produce reasonable agreement with prior forcing sets.Data sets that were instrumental in the EVA construction include the following: -The estimate of total SO 2 injection by Pinatubo of 18 ± 4 Tg (9 ± 2 Tg S) from the total ozone mapping spectrometer (TOMS) instrument (Guo et al., 2004).
-Aerosol extinction from the CCMI (Arfeuille et al., 2013).Here, we have used the CCMI data provided at http://www.pa.op.dlr.de/CCMI/CCMI_SimulationsForcings.html, specifically the files produced for use in the ECHAM6 model.We focus primarily on the EXT at 550 nm (EXT 550 ).AOD is computed based on the vertical integral of EXT above the climatological tropopause.Post-Pinatubo AOD and EXT anomalies are computed by subtracting an estimate of the background aerosol levels, based on the mean EXT field from the years 1999 and 2000.
-Estimates of aerosol effective radius are provided in the Sato/GISS and CU13 data sets.Approximate values of effective radius were also computed for the CCMI data set, by scaling the provided mean radius under the assumption of a log-normal distribution with an effective standard deviation of 1.2.
Input data, specifying the stratospheric sulfur injection properties for a number of volcanic eruptions, were collected from a range of sources and summarized in Table 1.The ice-core-derived sulfate aerosol loading estimates of Gao et al. (2008) were used to produce SO 2 injection estimates for the great Tambora eruption of 1815, as well as the eruptions of Agung (1963), Fuego (1974), and El Chichón (1982).For the 1991 eruption of Pinatubo, we use the estimated total SO 2 injection of 18 Tg (9 Tg S) from the TOMS instrument (Guo et al., 2004).Estimates of stratospheric SO 2 injection for a number of relatively smaller eruptions in the 2000s were taken from Brühl et al. (2015), based on the MIPAS SO 2 retrievals described by Höpfner et al. (2015).(Brühl et al., 2015) a Hemispheric asymmetry ratio defined as the estimated ratio of NH / SH AOD or stratospheric aerosol loading.b Taken as the ratio of AOD NH / AOD SH for the first 2 years after each eruption from the reconstruction of Stothers (2001).c Taken as the ratio of AOD NH / AOD SH for the first 2 years after the eruption from the CCMI reconstruction (Arfeuille et al., 2013;Thomason and Peter, 2006).

Global mean AOD and effective radius
Aerosol optical depth at λ = 550 nm is simulated by EVA making use of the assumption that it can be a simple function of the stratospheric sulfate mass: Stothers ( 1984) introduced a simple linear scaling between global sulfate aerosol mass and global mean AOD.Such scalings have a physical foundation (e.g., Charlson et al., 1992) and have long been used to study tropospheric aerosols.A linear relationship between sulfate mass and AOD has been used to convert the sulfate mass estimates of Gao et al. (2007) into radiative properties (Schmidt et al., 2011) and is implicit in the linear scaling of ice core sulfate flux to AOD used by CU13 for eruptions of Pinatubo magnitude and smaller.We apply this same assumption hereafter (up to a threshold M * ; see Sect.3.6) using a scaling factor A: Time evolution of sulfate mass is emulated in EVA using a chemical box model framework.Bluth et al. (1997) introduced a simple single-box model of stratospheric aerosol evolution.Changes in stratospheric sulfate aerosol are controlled by injections of SO 2 into the box, subsequent conversion of SO 2 into sulfate aerosols, and loss of sulfate aerosols to the troposphere.Describing the production and loss of sulfate mass (M SO 4 , all masses in Tg S) as a function of the SO 2 mass (M SO 2 ) and characteristic timescales τ prod and τ loss , respectively, the time tendency of This equation describing the time evolution of M SO 4 can be solved analytically, for example, for a single pulse injection of M SO 2 at time t 0 : Equations ( 2) and ( 4) can be used to emulate the observed global mean AOD evolution after the Pinatubo eruption.Using the best estimate of a total SO 2 injection of 9 Tg S (Guo et al., 2004), the parameters A, τ prod , and τ loss can be determined based on comparison with the CCMI global mean AOD anomaly time series (Fig. 1).Given the uncertainties in the CCMI AOD in the first months after the Pinatubo eruption due to the gap in the SAGE II observations due to saturation effects (Thomason and Peter, 2006), we have based our fit on the CCMI AOD beginning in July 1992, when SAGE II retrievals of the full tropical stratosphere resumed.Therefore, the fit is not strongly constrained by the peak of global mean AOD in the CCMI reconstruction, but rather by the shape of the AOD decay.Best fit is achieved with values of A = 0.0364, τ prod = 180 days, and τ loss = 330 days.The resultant value of A is comparable to that suggested by Stothers (1984), whose scaling factor is 0.0267 when converted into units of mass S rather than mass sulfate aerosol, and the stratospheric loss rate (τ loss = 330 days) is consistent with the decay rate of AOD after Pinatubo noted by other researchers (e.g., Bluth et al., 1997).The best fit production rate, τ prod = 180 days, is perhaps surprising, given its discrepancy from the observed timescale of SO 2 decay, which is around 30-35 days (Bluth et al., 1992;Read et al., 1993).In their single-box aerosol model, Bluth et al. (1997) noted the resulting lag between peak observed AOD and the peak in modeled SO 4 mass when using an SO 4 production timescale equal to the observed SO 2 decay rate, and proposed that the rate of SO 4 aerosol increase may be limited by other chemical steps other than the destruction of SO 2 .We contend that global mean AOD may further depend on the timescale of spatial spreading of the aerosol cloud, since the impact of aerosol contained within a horizontally contained vertical column will be diminished due to shielding effects.Whatever the mechanism responsible, τ prod should be interpreted as an effective production timescale, which incorporates not only the chemical conversion of SO 2 to SO 4 but other processes which damp the rise in AOD.An important caveat of this construction is that the peak loading of the simulated SO 4 mass is significantly less than what would result from a complete conversion of SO 2 to SO 4 prior to any loss.For this reason, the scaling factor A is larger than what would be deduced if the peak SO 4 loading is assumed to be equal to the SO 2 injection (in Tg S).Finally, we note that since estimates of the mass of SO 2 injected by Pinatubo (9 ± 2 Tg S, from Guo et al., 2004) have an uncertainty of about 25 %, the uncertainty in the scaling factor A is at least this large.The global mean effective radius is computed based on a simple scaling argument.For a given mass of sulfate M SO 4 , which is distributed equally among N aerosols of radius r, M SO 4 ∼ N r 3 . (5) If N is constant, then particle radius scales as the one-third power of M SO 4 .This relationship is similar to that used by CU13, who based r eff on the one-third power of AOD 550 .This follows if r is linearly related to r eff , which is the case, for instance, if the aerosol particles are distributed lognormally by size (with a given shape parameter).In EVA, we set and find a scaling constant R that produces best agreement in terms of the peak global mean r eff reached after the Pinatubo eruption in the r eff time series of the observation-based reconstructions.Like CU13, we also set a minimum r eff of 0.2 µm.With a fit value of R = 0.78, the resulting EVA r eff time series shows reasonable agreement in peak magnitude compared to the Sato/GISS, CCMI, and CU13 reconstructions with peak values between 0.5 and 0.6 µm (Fig. 1b).Basing r eff directly on the sulfate mass in this way does not allow the reproduction of some observed features of the r eff evolution apparent in the Sato/GISS r eff reconstruction, including the lag of its peak compared to that of AOD 550 .However, given the simplicity of this approach, and the uncertainties in r eff estimates retrieved from satellite sensors, the scaling methodology appears to produce satisfactory results, and even if the Sato/GISS evolution of r eff is more realistic, it remains to be demonstrated that this difference has any detectable influence on the climate.

Spatiotemporal structure
The stratosphere can be separated into regions with distinct dynamical regimes (Plumb, 1996(Plumb, , 2002)).One important distinction is that of the relatively undisturbed tropical pipe -where mean residual motion is predominantly upward -from the extratropics, where wave breaking leads to quasi-horizontal mixing, motivating the term "surf zone" (McIntyre and Palmer, 1983).The structure of stratospheric trace-gas species are well reproduced by the so-called "leaky pipe" model of the stratosphere, which differentiates between the different dynamical regimes of the tropics vs. extratropics (Ray et al., 2010).Following these studies, and motivated also by the clear separation of tropical vs. extratropical stratospheric aerosol maxima following the 1991 Pinatubo eruption (Trepte et al., 1994, and Fig. 2), EVA uses a simple three-box representation of the stratosphere, separating it into three regions -equatorial, Northern Hemisphere (NH) extratropical, and Southern Hemisphere (SH) extratropicaland describes the stratospheric aerosol distribution as the superposition of three zonally symmetric, global-scale aerosol plumes.
The aerosol properties of each plume are defined using a static characteristic spatial structure based on extinction from the CCMI reconstruction following the 1991 Pinatubo eruption.To reproduce the observed meridional structure of AOD, we choose for simplicity Gaussian functions, fit as a function of the area-conserving coordinate sin(ϕ), where ϕ denotes latitude.The observed evolution of zonal mean AOD 550 after Pinatubo (Fig. 2a) shows a clear separation between the three regions, with an initial growth of AOD 550 within the equatorial region and later peaks in the NH and SH.We base the structure of the equatorial shape function on the CCMI AOD 550 averaged over the first 3 months after Pinatubo, before much transport to the extratropics occurred, and fit a Gaussian function to the central portion of the AOD 550 (Fig. 2b).(While the magnitude of AOD 550 in the first months is assumed to be highly uncertain during the first months due to saturation effects, we necessarily assume here that the latitudinal structure of the CCMI AOD is reasonably realistic.)Subtracting the Gaussian fit of the equatorial plume from the 4-year mean AOD 550 isolates the spatial structure of the remaining extratropical regions (Fig. 2c).We fit the observed extratropical AOD structure with Gaussian functions centered at 45 • with a width of 14 • .At high latitudes, where the impact of our choice of the sin(ϕ) coordinate is strongest, the EVA fit shows some disagreement with the CCMI AOD, but much better agreement with the older Sato/GISS reconstruction (not shown).Given that SAGE observations are generally limited to ±60 • and high-latitude values are extrapolated, and that high-latitude aerosol often resides at heights below 15 km which can be difficult for satellite sensors to retrieve (Ridley et al., 2014), it does not appear warranted to change the sin(ϕ) fitting procedure adopted at lower latitudes.
The three horizontal shape functions are normalized by their respective global area-weighted mean.In this way, multiplication of a horizontal shape function by the sulfate mass in the region gives, for each latitude, a representation of the local vertically integrated sulfate mass density (kg km −2 ).
The time evolution of the AOD 550 for each region is based on expanding the single-box model of the stratosphere of Sect.3.2 into a three-box model.In addition to injections of SO 2 and loss of sulfate through cross-tropopause transport -as in the global single-box model -sulfate is transported between the three boxes, with time constants τ mix and τ res defining the rates of two-way mixing and poleward residual circulation, respectively.For example, the rate of change of sulfate mass in the NH region is related not just to SO 2to-SO 4 conversion and loss but also to mixing, which is related to the difference of mass between the NH and equatorial plumes, and the transport of mass from the equatorial plume due to the residual mass circulation of the Brewer-Dobson circulation: Similar expressions are used for the SH and equatorial regions.In the EVA module code, the differential equations describing the time tendency of SO 4 within each region are computed numerically through a forward Euler method: Each eruption is treated as an instantaneous injection of SO 2 into one of the three boxes, with the injection region based on the latitude of the volcano, and latitudinal boundaries set to ±25 • based on satellite-based estimates of the edges of the stratospheric tropical pipe (Neu et al., 2003).In this way, the impact of multiple eruptions occurring with overlapping time periods of impact can be easily handled, as each eruption simply adds to the pre-existing sulfur loading.
The observed AOD 550 after Pinatubo shows clear influence of the seasonal cycle of stratospheric transport, with www.geosci-model-dev.net/9/4049/2016/Geosci.Model Dev., 9, 4049-4070, 2016  maximum extratropical AOD found in the winter and spring seasons of each hemisphere, qualitatively consistent with the seasonal cycle of the Brewer-Dobson circulation which maximizes in the winter months (Holton et al., 1995).We therefore vary the τ mix and τ res parameters with calendar month m, using simple sinusoidal relationships, such as where the factor B describes the amplitude of the seasonal variation.Mixing and residual transport in the NH are thus strongest in January and weakest in July, and have an annual average equal to τ mix and τ res , respectively.Mixing and residual circulation in the SH are phase shifted by 6 months with maximum mixing in July and minimum in January.The mixing and residual transport timescales are equal for both hemispheres, and in the long term, produce relatively even partitioning of sulfate between the NH and SH.A method to reproduce hemispherically asymmetric aerosol evolution is described in Sect.3.6.
Fitting the parameters B, τ mix , and τ res was achieved by minimizing the root mean square residuals of the EVA AOD 550 with the CCMI AOD 550 field for Pinatubo.Owing again to the larger uncertainties in satellite-based aerosol properties in the initial months after Pinatubo due to saturation effects, the fitting procedure ignored the initial 12 months between the latitudes between 20 • S and 20 • N. Best fit was achieved with values of τ mix = 15 months, τ res = 17 months, and B = 0.75, resulting in good agreement with CCMI AOD 550 evolution (Fig. 3).The AOD 550 evolution of EVA lacks the fine detail of the CCMI data set, but reproduces the general spatiotemporal structure, including the double peak structure in the SH midlatitudes, and the gradual shift from strongest AOD 550 in the tropics to the extratropics with time.
In the vertical dimension, observations show that volcanic aerosol extinction values were strongly peaked in the tropics at about 22 km in the first months after the Pinatubo eruption, and thereafter spread somewhat with height, with the peak shifting slightly downwards (Arfeuille et al., 2013).Transport of aerosol from the equatorial region to midlatitudes and high latitudes was episodic in the first year after Pinatubo with contributions from both the lower and upper branches of the Brewer-Dobson circulation (Trepte et al., 1993).Horizontal transport of passive dynamical tracers in the midlatitude stratosphere takes place predominantly along isentropic surfaces, due to large-scale mixing processes resulting from wave activity.Aerosol transport is complicated by additional processes, including gravitational settling and anomalous vertical motions resulting from the heating of the local atmosphere by the absorption of longwave radiation by the aerosols themselves (Rogers et al., 1999).Nonetheless, the variation of peak aerosol extinction in the 4-year post-Pinatubo mean follows isentropic surfaces (derived from ERA-Interim reanalysis data; Dee et al., 2011) reasonably closely in the midlatitudes of the summer hemisphere, with the 430 K potential temperature surface corresponding well with the vertical peak of mean EXT (Fig. 4a).The correspondence between the aerosol loading and isentropic (potential temperature) surfaces breaks down in the high latitudes during winter, where diabatically driven vertical motion within the polar vortex leads to downwelling of air parcels with respect to potential temperature (Manney et al., 1994;Tegtmeier et al., 2008).
To reproduce the variation of the vertical peak of extinction with latitude in the extratropics, we define a vertical "center line" based on climatological potential temperature.While potential temperature does not perfectly follow the observed aerosol peak at high latitudes, by linking the vertical distribution of the aerosol to a temperature-based (rather than mass-or geopotential-based) vertical coordinate, it may be more suitable for application in much warmer or colder climates.Specifically, we define the center line at midlatitudes and high latitudes (ϕ > 45 • ) from the summer 430 K potential-temperature surface for each hemisphere (July for the NH, January for the SH), and the annual mean 430 K potential-temperature surface for −45 • < ϕ < 45 • .This empirically defined center line shows reasonable agreement with the observed vertical peak in aerosol extinction in the extratropics (Fig. 4a).In the tropics, the vertical position of the plume is given a vertical offset from the center line.
The vertical shape of extinction in the equatorial region is based on a Gaussian fit of the 4-year Pinatubo mean CCMI extinctions with a vertical width of 2.25 km and a vertical offset of 2.75 km (Fig. 4b).In the extratropics, a Gaussian fit is produced, centered on the defined center line with a width of 2.825 km (Fig. 4c).The vertical shape functions, defined internally on a 1 km vertical grid, are normalized by their vertical sum.Therefore, multiplication of the AOD at each latitude by the normalized vertical shape function produces a profile of aerosol extinction (per kilometer).
The resulting EVA aerosol extinction structure shows good agreement with the 4-year mean CCMI Pinatubo observations (Fig. 5).The use of Gaussian vertical and horizontal shape functions results in fairly good reproduction of the strong extinction gradient at the tropopause.
The spatiotemporal structure of aerosol effective radius is based on the simulated sulfate mass time series for each latitude (i.e., a function of the three-box model and the horizontal shape functions).Effective radius is computed for each latitude and time from Eq. 6, and assumed to be constant with height through the stratosphere.The resulting r eff field shows reasonable agreement with observation-based estimates (Fig. 6); however, it appears to produce peak values which are somewhat too large, and peak too early compared to the Sato/GISS reconstruction, as apparent also in the global mean (Fig. 1).

Wavelength-dependent optical properties
The wavelength-dependent optical properties EXT, SSA, and ASY are computable via Mie scattering theory given knowledge of the extinction at a specific wavelength, and the effective radius of an assumed log-normal size distribution.EVA utilizes look-up tables, computed from Mie theory, assuming a single-mode log-normal size distribution with width parameter σ = 1.2, which gives wavelength-dependent EXT scaling ratios (with respect to EXT at 550 nm), ASY, and SSA for varying effective radii, ranging from 0.2 to 1.3 µm, and for 29 wavelengths ranging from 0.2 to 100 µm.Therefore, given the EXT 550 and r eff at any point in space, EXT(λ), SSA(λ), and ASY(λ) are calculated from the lookup tables through bilinear interpolation.The use of σ = 1.2 is roughly consistent with that deduced from observations of the Pinatubo aerosol cloud, with Stenchikov et al. (1998) using σ = 1.25 and the CCMI reconstruction using σ = 1.2 in gap-filled regions when the SAGE measurements cannot be used to estimate sigma directly (Arfeuille et al., 2013).
Resulting zonal mean AOD at three sample wavelengths are shown in Fig. 7 and compared to the corresponding CCMI values.In general, EVA reproduces the decrease of AOD with increasing wavelength from the visible to the near infrared.SSA and ASY at four wavelengths, at 20 km, are shown in Figs. 8 and 9, and again reproduce reasonably well the variation of magnitude with changing wavelength.Some differences in the spatial structure of SSA are apparent between the EVA and CCMI reconstructions, although it remains to be shown how important such structure is in the overall climate impact of the volcanic forcing.to zero in the 2000-2010 decade has been shown to lead to non-negligible biases in climate model simulations (Solomon et al., 2011).The background stratospheric aerosol layer is thought to be a result of a combination of factors, including the episodic but ubiquitous influence of relatively minor volcanic eruptions with small stratospheric injections, and the slow and steady influx of aerosols and precursors from the troposphere into the stratosphere (Kremser et al., 2016).
A non-zero background stratospheric aerosol forcing is parameterized within EVA by specifying a constant SO 2 injection value, chosen to produce a best fit between the CCMI and EVA global mean AOD time series for the year 2000, where the observed global AOD reached a minimum (Fig. 10).This simple procedure results in a background SO 2 injection estimate of 0.2 Tg year −1 .For comparison, using the coupled aerosol Solar-Climate Ozone Links chemistry climate model (SOCOL-AER) and CCMI retrievals, Sheng et al. (2015) have inferred a total net sulfur mass flux into the stratosphere of 0.19 Tg year −1 .The structure of the observed aerosol extinction in the meridional height plane during the 2000 minimum is shown in Fig. 11a.Maximum aerosol extinctions in the CCMI data set are in the high-latitude lowermost stratosphere.To best reproduce this structure, the background SO 2 injections are split evenly between the two extratropical plumes.The resulting structure of EVA background AOD forcing does not reproduce the lower stratosphere forcing due to the static shape functions of the EVA construction based on the Pinatubo aerosol evolution.Nonetheless, this strategy of background AOD in EVA does somewhat reproduce the meridional structure of the background AOD (Fig. 11c, d).

Hemispheric asymmetry
Some tropical eruptions can lead to relatively even partitioning of aerosol between the NH and SH.For example (discounting the likely small and short-lived impact of the August 1991 Cerro Hudson eruption), the June 1991 eruption of Pinatubo produced NH and SH AOD maxima of similar magnitude (Fig. 3).On the other hand, tropical erup-  1994 1995 1996 1997 1998 1999 2000 2001 2002 2003  tions like Agung (1963) andEl Chichón (1982) are characterized by substantial hemispheric asymmetry in stratospheric aerosol loading.Asymmetries in the radiative forcing, if large enough, may have a detectable climate impact, for instance, through their influence on the latitude of the intertropical convergence zone and subsequent changes in tropical monsoon patterns (Haywood et al., 2013;Oman et al., 2006).
Some degree of hemispheric asymmetry may be expected due to the timing of an eruption with respect to the seasonal cycle of stratospheric transport.One expects tropical eruptions which occur in NH winter (when transport to the NH extratropics is strongest) to produce higher NH aerosol loadings than those in summer.This expectation is reproduced by aerosol general circulation models (Toohey et al., 2011).The seasonal variation of stratospheric transport parameterized within EVA leads to a small degree of hemispheric asymmetry in the resulting AOD patterns depending on season of eruption, with highest asymmetry produced for tropical eruptions in August (AOD NH / AOD SH = 1.18) and February (AOD NH / AOD SH = 0.847).
The parameterized seasonal stratospheric transport of EVA qualitatively reproduces the observed SH bias of the AOD from the 18 February 1963 Agung eruption (Fig. 12), albeit with much weaker magnitude.For instance, based on ground-based optical measurements, Stothers (2001) estimated that the SH aerosol loading after Agung was 8 times larger than that of the NH.Furthermore, the hemispheric asymmetry produced by the seasonal stratospheric transport of EVA for El Chichón, for an 4 April eruption date, is much different than the observed strong NH asymmetry (Fig. 12).The degree of hemispheric asymmetry for tropical eruptions may be related to the particular synoptic-scale meteorological conditions at the time and place of the initial SO 2 injection, and therefore impossible to reproduce using a climatological transport parameterization.
For these reasons, an anomalous asymmetry factor is included in EVA to allow one to impose an observed asym-metry from data.Hemispheric asymmetry can be specified in the input file, based on direct observational estimates or estimated from the ratio of Greenland to Antarctic ice core sulfate records.When this ratio differs from the asymmetry produced by the seasonal mixing parameterization of EVA, a correction can be applied which attenuates the rate of mixing and transport to one hemisphere for a period of 18 months after the eruption.This correction ensures that the AOD in the midlatitudes and high latitudes shows a hemispheric ratio www.geosci-model-dev.net/9/4049/2016/Geosci.Model Dev., 9, 4049-4070, 2016 equal to that of the input value, while retaining the nominal seasonality of stratospheric dynamics.
Hemispheric asymmetry corrections are applied here to the EVA reproductions of the Agung, Fuego, and El Chichón eruptions (Table 1).In order to enhance comparability with the CCMI data set, correction factors for Agung and Fuego are based on the AOD reconstruction of Stothers (2001), and defined as the ratio of the 2-year post-eruption average of hemispheric AOD.This method produces values of 0.19 for Agung and 1.0 for Fuego.For El Chichón, we calculate a correction factor directly from the CCMI AOD 550 , again as the ratio of the 2-year AOD average for each hemisphere, resulting in a value of 1.5.Applying hemispheric asymmetry corrections in the reconstruction of AOD 550 for Agung and El Chichón produces better agreement with the CCMI estimates (Fig. 12).It should be noted here that the CCMI representations for these eruptions are highly uncertain.For Agung, the CCMI data are based on 2-D model results scaled by sparse ground-based measurements, while for El Chichón, it is based on an amalgam of very sparse ground-and airplanebased lidar and satellite observations at only the high latitudes (Thomason and Peter, 2006).

Nonlinear scaling for large eruptions
While a simple linear relationship between sulfate mass and AOD (Sect.3.2) is consistent with available observations, its applicability to eruptions larger than Pinatubo is quite uncertain.Modeling studies imply that for large eruptions, the maximum AOD produced is not a linear function of the injected sulfur (Timmreck et al., 2010).CU13 argued that above some threshold, AOD should scale as the two-thirds power of sulfate aerosol mass rather than linearly.Such a relationship is consistent with the results of an aerosol general circulation model simulating a large range of tropical volcanic eruptions (Metzner et al., 2014).We sketch a simple explanation for this relationship as follows: for a total mass of sulfur (M) distributed among N particles with uniform radius r, the mass of each particle, M/N, is proportional to the volume of each particle (V ) and therefore r 3 : The AOD is proportional to the total cross-sectional area of the particles: Combining these equations, AOD can be written as a function of N and M: When new particles are formed in proportion to the mass of injected sulfur, N is proportional to M and the AOD scales linearly with M. On the other hand, if injected sulfur mass condenses onto pre-existing particles, N remains constant, and AOD scales with the two-thirds power of M. In EVA, we adopt a threshold-based implementation based on the approach of CU13.We retain this approach for very large eruptions, and scale sulfate to AOD based on two parameters, A 1 and A 2 , such that where the threshold M * is defined so as to make the relationship continuous: The two regimes -complete new particle formation from the injected SO 2 or complete condensation onto pre-existing particles -obviously represent the two extremes of possible sulfate aerosol evolution, and it seems likely that both processes take place in differing degrees for different eruptions.Other physical parameterizations, for instance, relating N to the logarithm of M are surely possible, and should be explored in future work.The present scheme retains consistency with the reconstruction of CU13, and has the advantage of simplicity, at least for the majority or eruptions for which AOD is a simple linear scaling of sulfate aerosol mass.Scaling considerations for extremely large eruptions should be 1960 1965 1970 1975 1980 1985 1990 1995 2000 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 1. Top: global mean AOD at 550 nm from the two reconstructions.Bottom two panels: zonal mean AOD at 550 nm as a function of latitude and time with log color scale.understood to be a major source of uncertainty in any volcanic forcing reconstruction.
In EVA, the A 1 term is based on the peak global mean AOD from the CCMI data set and the best estimate of the 9 Tg S injection from Pinatubo (Sect.3.2).We choose to define A 2 so as to best reproduce the aerosol forcing of CU13 for large eruptions, such as Tambora and Samalas.
Figure 13 shows the relationship between peak global mean AOD and estimated injected SO 2 for Pinatubo, Tambora, and Samalas.The linear scaling of CU13, which was based on AOD estimates from Sato/GISS and ice-core-based estimates of SO 2 injection, leads to a rather steep curve, which, when extrapolated to the estimated sulfur injection of Tambora, would have produced an estimate of global mean AOD of about 0.7.A much smaller AOD for Tambora (i.e., 0.45) and Samalas was produced by CU13 by incorporating the two-thirds power law, which they applied to eruptions larger than Pinatubo, producing the AOD values shown in Fig. 13.
The linear scaling used in EVA is significantly less steep than that of CU13, a result of the lower peak global mean AOD estimate for Pinatubo from CCMI compared to Sato/GISS, and the larger estimate of SO 2 injection from satellite sensors compared to the ice-core-derived estimate of CU13.Extrapolation of the linear scaling of EVA to larger injection magnitudes reproduces the CU13 estimates of AOD and r eff for Tambora rather well.A two-thirds power-law relationship is nonetheless implemented in EVA, which applies then to eruptions greater in magnitude than Tambora, and therefore reduces the impact of the Samalas eruption compared to using the linear relationship.

Modern era
AOD time series produced by EVA using the eruption history of Table 1 are shown in Fig. 14, and compared to the CCMI data set.The magnitude of global mean AOD of the major eruptions is well reproduced by EVA, with slightly larger AOD produced for Pinatubo (as discussed in Sect.3.2) and slightly lower AOD for Agung compared to the CCMI data, which relies on ground-based optical data to scale the AER model results for Agung and Fuego.A number of relatively minor eruptions before and after the El Chichón eruption (Bluth et al., 1997) -not accounted for in the eruption history used here -likely account for the complex temporal and spatial AOD evolution seen in the CCMI database in the years 1980-1990.Including the minor eruptions of 2000-2014 in the EVA input file improves the comparison over these years.This time period is examined in more detail in Fig. 15 1.Top: global mean AOD at 550 nm from the two reconstructions.Bottom two panels: zonal mean AOD at 550 nm as a function of latitude and time.Sarychev (2009) exceeds that seen in the CCMI data, and appears to persist longer.This is likely a result of the use of a single aerosol loss rate for all eruptions in EVA, based on the observed decay of the aerosol from Pinatubo.It is likely that the processes and related timescales are different for relatively smaller eruptions in the extratropics, but there is presently little understanding of the relationship between eruption strength and resulting aerosol lifetime and there is no consideration of this potential effect in EVA.Nonetheless, we note that the agreement between EVA and CCMI is decent at ∼ 60 • N, the northernmost limit of the CALIPSO measurements underlying the CCMI data, implying that the stronger apparent disagreement at polar latitudes may be due partly to the gap-filling procedure used in the construction of the CCMI data set.The agreement between EVA and CCMI over this time period is also largely dependent on the accuracy of the stratospheric injection estimates used.The estimates used here, based on MIPAS satellite SO 2 measurements, carry an uncertainty on the order of 20 % (Höpfner et al., 2015), with additional (unquantified) uncertainty arising from the separation of the stratospheric component of the total atmospheric injection (Brühl et al., 2015).If desired, agreement between the EVA-based results with the satellite records over the 2004-2014 period could be improved by careful adjustment of the input sulfur injections.

Tambora
Global mean AOD 550 and r eff produced by EVA for the 1815 eruption of Tambora are compared to the reconstruction of CU13 in Fig. 16.Using the linear scaling, based on the estimated peak AOD and sulfur injection of Pinatubo, the peak AOD of an estimated 55 Tg SO 2 injection by Tambora leads to a smaller global mean AOD in the EVA reconstruction compared to CU13.The EVA AOD peaks at approximately 0.35, while the CU13 reconstruction estimates a global mean AOD peak of ∼ 0.4.
It should be noted that with the Pinatubo-based linear relationship used in EVA, if a nonlinear, two-thirds power-law relationship were implemented with the same threshold used by CU13, i.e., applying to eruptions just greater in magnitude than Pinatubo, then the resulting AOD for Tambora would be significantly smaller than the present EVA estimate.
The aerosol effective radius produced by EVA for Tambora is similar, but slightly larger than the CU13 estimate.This difference comes about because in EVA, r eff is scaled according to the SO 4 mass, while in CU13, the AOD is used.Since in CU13, the AOD for Tambora is related to the twothirds power of mass, the muting of the AOD affects also the r eff , while in EVA, r eff grows linearly with SO 4 mass.There are unfortunately no observational estimates of the aerosol size distributions for any eruptions larger in magnitude than Pinatubo (1991), and therefore there is considerable uncertainty regarding the validity of the simple relationships used here and by CU13 for such large eruptions.Simulations of the Tambora eruption with an interactive aerosol model produce similar effective radii as those estimated by EVA (∼ 0.73 µm; Stoffel et al., 2015); future work with such models may help to constrain parameterizations of effective radius.
The zonal mean AOD 550 for Tambora from CU13 and from this work are shown in Fig. 17.While the magnitude of peak AOD and its temporal decay are similar in the two reconstructions, the EVA reconstruction provides a more realistic latitudinal distribution of AOD compared to the 4-band structure of CU13.The slight SH bias of the CU13 Tambora AOD 550 , inferred by CU13 from ice core records, is produced by EVA as a result of the parameterized seasonal transport (Sect.3.3) with no additional hemispheric correction.

Conclusions and outlook
EVA has a number of strengths, drawing on the advantageous aspects of previous volcanic aerosol reconstructions.Like the Stenchikov and CCMI data sets, EVA provides full field optical properties (EXT, SSA, ASY) as a function of wavelength, height, and latitude, allowing consistent implementation within different climate models.Like the Amman et al. (2003) and Gao et al. (2008) reconstructions, EVA provides a consistent treatment of all eruptions, avoiding gaps and discontinuities that are unavoidable in reconstructions based only on observations.EVA produces good agreement with satellite-based observations of the Pinatubo eruption, providing confidence in its ability to produce reasonably realistic forcing structure.For the 1815 Tambora eruption, the AOD produced by EVA is similar to prior reconstructions (see Sect. 4.2), and approximately in the center of the broad range of estimates produced by interactive stratospheric aerosol models (Zanchettin et al., 2016).With standard parameter settings, EVA therefore provides a middle-ofthe-road forcing estimate for given eruption magnitudes.At the same time, the construction of EVA allows for great flexibility.Different EVA reconstructions can be constructed using different histories of stratospheric sulfur injections from observations (Brühl et al., 2015;Carn et al., 2016;Höpfner et al., 2015;Neely and Schmidt, 2016) or from ice cores (e.g., Gao et al., 2008), thereby translating uncertainties in the volcanic emission estimates into aerosol forcing parameters.Adjusting EVA reconstructions to be consistent with updated satellite retrievals or with model results is achievable through modification of the parameter settings.This flexibility also makes EVA well suited for idealized studies: for instance, the impact of different uncertainties in volcanic aerosol properties could be tested through producing an ensemble of EVA forcing sets with an ensemble of parameter settings.The simplicity of EVA also ensures that it is fast and can provide forcing reconstruction instantaneously for any given past or future eruption.
By design, EVA is simple, and cannot reproduce all the features of aerosol forcing which are seen in observationbased reconstructions or aerosol general circulation model results.The three-box representation of stratospheric transport neglects the impact of the polar vortex, which creates a mixing barrier in the winter hemisphere and likely enhances aerosol loss in a seasonal manner.EVA also does not presently consider the height of stratospheric sulfur injection, which likely plays an important role in the timescale of cross-tropopause transport (Bluth et al., 1997).Similarly, EVA does not account for vertical variations in stratospheric dynamics, specifically the shallow and deep branches of the BDC, which have different strengths and seasonal variability.
For most purposes, inaccuracies in forcing due to the simple approach of EVA are likely small compared to uncertainties in our knowledge of the properties of past volcanic eruptions inferred from proxies like ice cores.Instead of attempting to perfectly reproduce observed aerosol properties, EVA makes it possible to pose the scientific question as to what aspects of the volcanic aerosol can produce a detectable climate response, thereby providing a means for deepening our understanding of the interaction of the stratospheric aerosol and climate.Future updates to EVA are planned, but in keeping with its original motivation, no attempt will be made to represent all aspects of the observations, but only those that can be demonstrated to have a detectable influence on the climate response.Updates will be motivated by new and updated observations and, to the extent they can reliably constrain remaining uncertainties, information from more complex aerosol models.In the meantime, in its present form, EVA is useful for a variety of purposes, including producing forcing for paleo-modeling simulations, e.g., within the Paleo-Modelling Intercomparison Project (PMIP) (Kageyama et al., 2016), and for idealized volcanic forcing experiments such as those within the Model Intercomparison Project on the climate response to volcanic forcing (VolMIP) (Zanchettin et al., 2016).Additional potential uses of EVA include decadal prediction simulations in the case of a major eruption (Timmreck et al., 2016), filling gaps in satellitebased forcing reconstructions, or in experiments aiming to assess what aspects of the stratospheric aerosol forcing can lead to a detectable climate response.

Code availability
EVA version 1.0 code, a user's manual, sample input data files, and driver scripts are included as a Supplement.Future updates and development will be coordinated through GitHub; see https://github.com/matthew2e/easy-volcanic-aerosol.

Figure 2 .
Figure 2. Definition of latitudinal shape functions.(a) The zonal mean CCMI AOD anomaly at 550 nm as a function of latitude and time for 4 years.(b) CCMI AOD anomaly averages over the 3 months (gray) after the Pinatubo eruption, normalized by its maximum value.A Gaussian fit to the equatorial portion of the 3-month average (blue) is used to define the latitudinal shape function for the equatorial plume of EVA.(c) The residual of the CCMI 4-year mean AOD minus the EVA equatorial shape function (gray) is used to define the extratropical EVA shape functions (blue).

Figure 4 .
Figure 4. Definition of EVA vertical shape functions.(a) CCMI aerosol extinction averaged over the first 4 years after the Pinatubo eruption, as a function of latitude and height.The climatological tropopause is shown in black and climatological July potential temperature (K) surfaces shown in gray for values as labeled.The EVA vertical center line (z center ), as defined in text, is shown in blue.(b) Normalized CCMI 4-year post-Pinatubo average extinction profiles as a function of altitude with respect to the vertical center line for the tropical (−15 to 15 • N) latitudes (gray), and the Gaussian fit defining the EVA vertical shape function for the equatorial plume (blue).Panel (c) is the same as (b) for extratropical (30-60 • ) profiles.

Figure 7 .
Figure 7. AOD after the Pinatubo eruption for selected wavelengths from (left) CCMI and (right) EVA.

3. 5 Figure 8 .
Figure 8. SSA after the Pinatubo eruption at 20 km for selected wavelengths from (left) CCMI and (right) EVA.Note different color scale for the lowermost panels.

Figure 9 .
Figure 9. ASY after the Pinatubo eruption at 20 km for selected wavelengths from (left) CCMI and (right) EVA.

Figure 10 .
Figure10.Global mean AOD for the 1994-2004 period from the EVA single-box approach and the CCMI, Sato/GISS, and CU13 reconstructions.

Figure 11 .Figure 12 .
Figure 11.Top: aerosol extinction at 550 nm averaged over the year 2000 from (left) the CCMI observation-based reconstruction and (right) EVA.Bottom: AOD at 550 nm through the year 2000 from (left) CCMI and (right) EVA.

Figure 13 .
Figure 13.Relationships between eruptive stratospheric sulfur injections and peak global mean AOD 550 from the CU13 reconstruction and EVA. Green crosses show the peak AOD values and injection estimates values from the CU13 reconstruction for Pinatubo, Tambora, and Samalas.The green line shows the linear relationship which would be deduced from the Pinatubo data used in CU13, extrapolated to all injection magnitudes.The blue cross shows the peak AOD and sulfur injection estimates from satellite sensors used in EVA to construct a linear relationship, shown by the blue solid line.The blue dashed line shows the two-thirds power-law relationship used in EVA for eruptions larger in magnitude than Tambora.

Figure 14 .
Figure14.Aerosol optical depth over the period 1960-2015 from the CCMI data set and from EVA, using the volcanic eruption history of Table1.Top: global mean AOD at 550 nm from the two reconstructions.Bottom two panels: zonal mean AOD at 550 nm as a function of latitude and time with log color scale.

Figure 15 .
Figure15.Aerosol optical depth over the period 2004-2012 from the CCMI data set and from EVA, using the volcanic eruption history of Table1.Top: global mean AOD at 550 nm from the two reconstructions.Bottom two panels: zonal mean AOD at 550 nm as a function of latitude and time.

Figure 16 .Figure 17 .
Figure 16.Global mean AOD at 550 nm and aerosol effective radius (r eff ) for the 1815 Tambora eruption from the CU13 reconstruction and EVA.

Table 1 .
Stratospheric sulfur injection estimates used in this work.