A perturbed parameter model ensemble to investigate Mt . Pinatubo ’ s 1991 initial sulfur mass emission

We have performed more than 300 atmospheric simulations of the 1991 Pinatubo eruption using the AER 2D sulfate aerosol model to optimize the initial sulfur mass injection as a function of altitude, which in previous modeling studies has often been chosen in an ad hoc manner (e.g., by applying a rectangular-shaped emission profile). Our simulations are generated by varying a four-parameter vertical mass distribution, which is determined by a total injection mass and a skew-normal distribution function. Our results suggest that (a) the initial mass loading of the Pinatubo eruption is approximately 14 Mt of SO2; (b) the injection vertical distribution is strongly skewed towards the lower stratosphere, leading to a peak mass sulfur injection at 18–21 km; (c) the injection magnitude and height affect early southward transport of the volcanic clouds as observed by SAGE II.


Introduction
The eruption of Mt.Pinatubo on 15 June 1991 injected large amounts of sulfur dioxide into the stratosphere.It perturbed the radiative, dynamical and chemical processes in the Earth's atmosphere (McCormick et al., 1995) and caused a global surface cooling of approximately 0.5 K (Dutton and Christy, 1992).The Pinatubo eruption serves as a useful analogue for geoengineering via injection of sulfur-containing gases into the stratosphere (Crutzen, 2006;Robock et al., 2013).Therefore, modeling volcanic eruptions advances our knowledge not only of the eruptions themselves on weather and climate, but also potential impacts of stratospheric sulfate geoengineering.
The uncertainties in determining the initial total mass and altitude distribution of SO 2 released by Pinatubo remain high.Stowe et al. (1992) deduced a mass of 13.6 megatons of SO 2 based on the aerosol optical thickness observed by the Advanced Very High Resolution Radiometer (AVHRR).By analyzing SO 2 absorption measurements from the Total Ozone Mapping Spectrometer (TOMS) satellite instrument, Bluth et al. (1992) estimated an initial mass loading of approximately 20 Mt of SO 2 .This study was later reevaluated by Krueger et al. (1995), who determined a range of 14-28 Mt emitted by Pinatubo, given the large retrieval uncertainties associated with TOMS.Later, Guo et al. (2004) constrained this range to 14-22 Mt of SO 2 .Besides the total emitted mass, the altitude distribution of the SO 2 emission is also not well constrained.The only available measurements with vertical resolution of SO 2 in the stratosphere during the Pinatubo period have been made by the Microwave Limb Sounder (MLS) in September 1991 (Read et al., 1993), which unfortunately only started its mission 3 months after the eruption.Given the lack of measurements in the period immediately following the Pinatubo eruption, modeling studies of Pinatubo (e.g., Weisenstein et al., 1997;Timmreck et al., 1999;SPARC, 2006;Heckendorn et al., 2009;Niemeier et al., 2009;Toohey et al., 2011;Aquila et al., 2012;English et al., 2013;Dhomse et al., 2014) have employed very different mass loadings, emission altitudes and vertical mass distributions, which leads to biases in the local heating and consequently in the dynamical response and time evolution of the stratospheric aerosol burden.These uncertainties, in addition to model-specific artifacts, make it difficult to accurately simulate the Pinatubo eruption.
Here, we attempt to provide a solution to the problems outlined above.We use the AER 2-D size-bin resolving (also called sectional or spectral) sulfate aerosol model (Weisenstein et al., 1997), which participated in an international aerosol assessment (SPARC, 2006), and was one of the bestperforming stratospheric aerosol models (in terms of comparing SO 2 , aerosol size distributions and extinctions with observations) under both background and volcanic conditions.We present results from more than 300 atmospheric simulations of the Pinatubo eruption based on different combinations of four emission parameters, namely the total SO 2 mass and a three-parameter skew-normal distribution of SO 2 as a function of altitude.We calculate aerosol extinctions from all of the simulations and compare them with Stratospheric Aerosol and Gas Experiment II (SAGE II) measurements (Thomason et al., 1997(Thomason et al., , 2008)).Such a head-on approach is currently impossible for global 3-D models due to computational expenses.The purpose of this work is to provide a universal emission scenario for global 3-D model simulations.To this end, we optimize the emission parameters such that the resulting SO 2 plume, aerosol size distributions, aerosol burdens and extinctions match balloon-borne, satellite and lidar measurements.We repeat two simulations using the 3-D SOCOL-AER aerosol-chemistry-climate model (Sheng et al., 2015) as a consistency check in a more complex model.In Sect. 2 we describe the model and the experimental design of our Pinatubo simulations.Section 3 compares the Pinatubo simulations with the observations, and conclusions follow in Sect. 4.

AER 2-D sulfate aerosol model
The AER 2-D sulfate aerosol model participated in an international aerosol assessment (SPARC, 2006), in which it was compared with satellite, ground lidar and balloon measurements, as well as with other 2-D and 3-D aerosol models, and subsequently recognized as one of the best existing stratospheric aerosol models with respect to SO 2 , aerosol size distributions and extinctions under both background and volcanic conditions.The model represents sulfuric acid aerosols (H 2 SO 4 /H 2 O) on the global domain from the surface to about 60 km with approximately 9.5 • horizontal and 1.2 km vertical resolution.The model is driven by yearby-year wind fields and temperature from Fleming et al. (1999), which were derived from observed ozone, water vapor, zonal wind, temperature, planetary waves, and quasi-biennial oscillation (QBO).The model chemistry includes the sulfate precursor gases carbonyl sulfide (OCS), sulfur dioxide (SO 2 ), sulfur trioxide (SO 3 ), sulfuric acid (H 2 SO 4 ), dimethyl sulfide (DMS), carbon disulfide (CS 2 ), hydrogen sulfide (H 2 S) and methyl sulfonic acid (MSA).The model uses pre-calculated values of OH and other oxidants from Notholt et al. (2005).Photodissociation and chemical reactions are listed in Weisenstein et al. (1997) and their rates are updated to Sander et al. (2011).The particle distribution is resolved by 40 size bins spanning wet radii from 0.39 nm to 3.2 µm by volume doubling.Such a sectional approach was proven to be more accurate in representing aerosol mass/extinctions compared to prescribed unimodal or multimodal lognormal distributions (Weisenstein et al., 2007).The sulfuric acid aerosols are treated as liquid binary solution droplets.Their exact composition is directly derived from the surrounding temperature and humidity according to Tabazadeh et al. (1997).Microphysical processes in the model include homogeneous nucleation, condensation/evaporation, coagulation, sedimentation, as well as tropospheric rainout/washout.These processes determine the evolution of the aerosol concentration in each size bin and thereby the entire particle size distribution.Operator splitting methods are used in the model with a time step of 1 hour for transport, chemistry, and microphysics, and 3 min sub-steps for the microphysical processes that exchange gasphase H 2 SO 4 with condensed phase, and 15 min sub-steps for the coagulation process.For more detailed descriptions of chemistry and microphysics in the model we refer to Weisenstein et al. (1997Weisenstein et al. ( , 2007)).

Coupled 3-D aerosol-chemistry-climate model
We employ the coupled aerosol-chemistry-climate model SOCOL-AER (Sheng et al., 2015) in order to verify the consistency between a 2-D model forced with observed dynamics and a 3-D free-running model.SOCOL-AER couples the size resolved AER 2-D microphysical model into the chemistry-climate model SOCOL (Stenke et al., 2013) with interactive aerosol radiative forcing.In this study we use the T31 horizontal resolution (3.75 • ×3.75 • ) and 39 vertical levels (from surface to 0.01 hPa) with nudged quasibiennial oscillation.Transport is calculated every 15 min, whereas chemistry, microphysics and radiation are calculated every 2 hours with 40 sub-steps (3 min) for the microphysics.This model has been well validated by comparing calculations with sulfur-containing gases, aerosol extinctions at different wavelength channels (from 525 nm to 5.26 µm), and aerosol size distributions from satellite and in situ observations.It has been used to study the global atmospheric sulfur budget under volcanically quiescent conditions and moderate volcanic eruptions such as the 2011 Nabro eruption.A detailed description of SOCOL-AER is presented by Sheng et al. (2015).

Experiments
We have simulated the Pinatubo-like eruption by injecting SO 2 directly into the stratosphere.In the 2-D model, the injection is immediately mixed zonally, and takes place in the latitude band 5 • S-14 • N, which is an approximation to the observed rapid zonal transport of the SO 2 cloud derived from satellite measurements (Bluth et al., 1992;Guo et al., 2004).The lack of zonal resolution is clearly a deficiency of our approach, but since SO 2 removal/conversion rate (efolding time) is sufficiently slow (τ ∼ 25 days) and the zonal transport around the globe sufficiently fast (τ ∼ 20 days) (Guo et al., 2004), a zonal-mean description is a reasonable approximation.Also, the spaceborne aerosol data are typically provided as zonal averages.We examined three cases of total mass, namely 14, 17 and 20 Mt of SO 2 .The injection height extends from near the tropical tropopause (17 km) to 30 km.The vertical mass distribution is then represented by M tot F (z) where M tot is the SO 2 mass magnitude in units of megaton (Mt) and ) is a vertical distribution function of altitude z ∈ [17 km, 30 km] with a skew-normal distribution f (z) given by (Azzalini, 2005) Figure 1 shows a few examples of F (z).The location parameter µ depends on available model levels and determines the altitude where the maximum of the emitted SO 2 cloud is located when there is no skewness.The skewness or asymmetry of the curve increases when |α| increases and vanishes when α = 0 (normal distribution).A negative α drives the location of the maximum SO 2 emission to lower altitudes, while a positive α to higher altitudes.The scale parameter σ indicates how much dispersion takes place near the maximum; that is, it determines the width or standard deviation of the asymmetric bell-shaped curve.
The four parameters M tot , µ, σ and α enable representation of a substantial space of SO 2 distributions, whose evolution is computed forward in time (taking into account the transport and comprehensive chemical and microphysical processes), in order to compare with satellite and in situ data.We simulate the following cases in detail: which results in 324 different scenarios.The choice of the boundaries for this set of scenarios is already based on exploratory simulations.For example, based on the results of our 2-D model, it does not make sense to consider total masses M tot > 20 Mt, since no choice of the other three parameters would allow to reconcile the model results with the observations.Similarly, skewness α > 0 leads to more biased model results, because the skew towards higher altitudes cannot be offset by lower M tot .In addition to the above 324 simulations, we consider another two scenarios, which are adopted in modeling studies of Pinatubo: (1) Box14Mt has a uniform ("Box") profile, which is similar to Dhomse et al. (2014) and the simulation "CONTROL_HIGH" in Aquila et al. (2012), injecting the SO 2 mass homogeneously along altitudes (shown in Fig. 1); (2) SPARC20Mt is the reproduction of the Pinatubo simulation conducted in SPARC ( 2006), which injects 20 Mt of SO 2 and has a vertical profile "SPARC" shown in Fig. 1.A selected list from the 326 2-D simulations is summarized in Table 1, in which the specific choice of the four parameters for each scenario is provided.The score and ranking of these scenarios are discussed later in the text.
Given the limitation of the 2-D approach, we further perform two 3-D Pinatubo-like simulations (R001 3-D and R153 3-D at the bottom of Table 1) using the coupled aerosolchemistry-climate model SOCOL-AER Sheng et al. (2015) to check the consistency between 2-D and 3-D approaches.Note that the location parameters used in the 3-D runs differ slightly from the corresponding 2-D runs (i.e., R001 and R153) due to different vertical model levels between the two models.

Metrics and data sets
To determine an optimal set of the emission parameters, we define four metrics (ScoreSO 2 , ScoreBurden, ScoreOPC and ScoreExt) based on these four measurements sets described above, and rank all of our 326 simulations by a weighted score (ScoreWt) of the four metrics (see Table 1).ScoreSO 2 is calculated as the relative l 2 -norm (Euclidean norm) error with respect to the MLS measurements: where X is a 1-D vector of SO 2 mixing ratio in altitude (21, 26, 31, 36 and 41 km).The negative values of the MLS measurements are set to zero in the calculation.ScoreBurden is the average of the relative l 2 -norm errors with respect to HIRS (July-December 1991) and SAGE-4λ (January 1992-December 1993): where B t 1 is a 1-D (in time) vector of the aerosol burden for July-December 1991 and B t 2 for January 1992-December 1993.ScoreOPC -We first calculate the relative l 2 -norm errors with respect to the OPC measurements: where N is a 1-D vector of the cumulative particle number concentration in altitude (15-30 km).We then evaluate a quadratic mean (RMS): rmsOPC = RMS{errOPC r }, where r denotes four particle size channels (r > 0.01 µm, r > 0.15 µm, r > 0.25 µm and r > 0.5 µm).Finally, Score-OPC is obtained by averaging rmsOPC from October 1991 to December 1992.
ScoreExt -The uncertainty of SAGE is generally better than ∼ 20 % for 525 nm and ∼ 10 % for 1020 nm (see Fig. 4.1 in SPARC, 2006).Therefore, ScoreExt is weighted as one-third for 525 nm (ScoreExt525nm) and two thirds for 1020 nm (ScoreExt1020nm).We use the SAGE II observations between 18 and 30 km.The calculations for Score-Ext525nm and ScoreExt1020nm are similar to those in ScoreOPC.Latitude bands (50-40 • S, 30-20 • S, 5 • S-5 • N, 20-30 • N and 40-50 • N) take the place of the particle size channels.The temporal average is from January 1992 to December 1993.
Note that extinction coefficients in the lower stratosphere (18-23 km) have a much larger weight than those above 23 km because extinctions at 525 nm and 1020 nm at 18-23 km after the Pinatubo eruption (see Fig. 5) are one to several orders of magnitude larger than those above 23 km.We calculate the score by the relative Euclidean norm, therefore the scores above 23 km have a relatively small weight.
The overall score ScoreWt is weighted as follows: 16.7 % of the SO 2 score (ScoreSO 2 ), 16.7 % of the OPC score (ScoreOPC), 33.3 % of the global burden score (ScoreBurden), and 33.3 % of the aerosol extinction score (ScoreExt).The choice of the weighting is discussed below.
MLS detected residual SO 2 in the stratosphere approximately 100 days after the eruption.The uncertainty of ScoreSO 2 is likely larger than ScoreBurden and ScoreExt due to uncertain OH fields.An assumed uncertainty in OH fields of 10 % (e.g., Prinn et al., 2005) translates into an uncertainty of 30 % in SO 2 at ∼ 90 days after the eruption.Moreover, ScoreOPC also has less weight than ScoreBurden and ScoreExt because of the small temporal and spatial sample size of the balloon-borne OPC measurements, which are not conducted very frequently (a maximum of two measurements per month after the Pinatubo eruption) and located only above Laramie.
ScoreBurden uses the HIRS-derived data up to December 1991 and the SAGE-derived data afterwards.During the first 6 months after the Pinatubo eruption, the SAGE II instrument was largely saturated in the tropical region (Russell et al., 1996;Thomason et al., 1997;SPARC, 2006;Arfeuille et al., 2013), and therefore the aerosol mass retrieved from SAGE II during this period very likely underestimates the initial loading significantly.The SAGE-4λ data set corrects for this deficiency by filling observational gaps by means of Lidar data.However, Lidar-derived extinctions are generally lower than SAGE II below 21 km (SPARC, 2006), and are not located in the equatorial region (see Fig. 3.7 in SPARC, 2006), where maximum mass loadings are expected.Therefore, SAGE II gap-filled data probably remain as a lower limit after the eruption.Conversely, HIRS measurements represent an upper limit since they account for the entire aerosol column including the troposphere.This may explain the considerable difference between SAGE II and HIRS during the first 6 months after Pinatubo (see Fig. 3).After this period, the aerosol mass in the extratropics contributes more to the global value than that in the tropics because the volcanic cloud starts to spread out from the tropics in November 1991 (see Fig. 5 of Baran and Foot, 1994).HIRS loses its sensitivity at mid/high latitudes where there is a contribution from errors in the background signal (Baran and Foot, 1994).As shown in Fig. 3, a visible increase of the HIRS-derived global burden begins after December 1991, and the noises in HIRS become more pronounced after March 1992.On the other hand, SAGE II, as an occultation instrument, becomes more reliable when the stratosphere starts to be sufficiently transparent after December 1991, particularly in mid latitudes.Therefore, ScoreBurden uses the HIRS-derived data up to December 1991 and the SAGE-derived data afterwards, with an overall uncertainty of 20 %.ScoreExt uses the SAGE II measurements from January 1992 to exclude the most saturated phase of SAGE II.

Scoring table
Table 1 shows the scores of selected scenarios, sorted according to the weighted rank ("RankWt" in the next to last column).The rank computed by the arithmetic average of the four scores is also provided ("RankAvg" in the third column from the right).The top 20 scenarios reveal that the total injection mass (M tot ) is 14 Mt of SO 2 , 70-80 % of which is below 24 km, and its maximum is likely between 18-21 km with 3-4 km width (scale parameter σ ).Location parame-ters µ larger than 21 km are skewed towards a lower altitude (negative α).These sort of vertical profiles provide a range for the parameters of the optimal vertical distribution: µ = 20.66 ± 1.79 km, σ = 3.33 ± 0.72 km and α = −0.8∓ 0.77.Two examples (scenarios R001 and R010 marked in Table 1) are shown in Fig. 1.The ranking based on "RankAvg" differs slightly from "RankWt"; however the set of the best scenarios found in "RankAvg" is consistent with "RankWt" despite the distinct weighting schemes.The worst scenarios ("RankWt" ≥ 322) in Table 1 are those with 20 Mt SO 2 injection mass and highest location parameters (µ = 29.55 km).The scenarios such as Box14Mt and R153 rank much more poorly than the optimal scenarios, although their injection mass is the same, because their vertical profiles (shown in Fig. 1) inject over 50 % mass above 23-24 km.The scenario R033 has the same vertical profile as R001, but injects 17 Mt SO 2 .SPARC20Mt emits 20 Mt SO 2 and ranks at 214 in Table 1, although its vertical profile is close to the optimal scenarios (about 10-20 % more mass above 23 km).This implies that emitting above 17 Mt SO 2 is very likely an overestimation.
The optimal vertical profiles found in Table 1 are generally consistent with the earlier volcanic plume studies of Fero et al. (2009) and Herzog and Graf (2010).Fero et al. (2009) showed that the SO 2 plume from the 1991 Pinatubo eruption originated at an altitude of ∼ 25 km near the source and descended to an altitude of ∼ 22 km as the plume moved across the Indian Ocean.Herzog and Graf (2010) suggested that initially SO 2 from a co-ignimbrite eruption (such as Pinatubo) that was forced over a large area, may reach above 30 km but the majority of SO 2 would then collapse or sink back to its neutral buoyancy height (15-22 km) (see Fig. 1 in their paper).
We discuss in detail nine scenarios (R001, R010, R017, R033, R153, Box14Mt, SPARC20Mt, R001 3-D and R153 3-D).R001 represents the overall optimal scenario.R010 ranks first in the ScoreExt and third in the ScoreBurden, as an example of scenarios with high rankings in the extinction and aerosol burden scores.R017 matches best the OPC measurement, but has poorer scores in the other criteria than R001 and R010.R086 has a vertical profile similar to R001 (see Fig. 1), and agrees the best with the SO 2 observations.However, this scenario fails to match other observations due to its abundant initial injection of 20 Mt SO 2 .R033 emitted 17 Mt SO 2 with the same vertical profile of R001, and ranks third in the ScoreSO 2 but poorly among other scores, which shows a performance similar to R086.Here we will focus on R033 for later discussion.R153 and Box14Mt (with RankWt 94) inject the same sulfur mass as in R001, but use different vertical profiles (maximum injection mass of R153 is located at ∼ 26 km).SPARC20Mt turns out to be a bad representation, which reproduces the previous simulation conducted in SPARC (2006).The two 3-D scenarios R001 3-D and R153 3-D correspond to the 2-D scenarios R001 and R153, respec-   (Read et al., 1993).
tively.The scores of the 3-D runs are similar to the corresponding 2-D ones.

Matching SO 2
Figure 2 compares the modeled SO 2 with MLS measurements in September 1991.The scenario R001 captures the measured SO 2 profile, and only underestimates the measured maximum SO 2 mixing ratio near 26 km by about 20 %.SO 2 modeled by R033 agrees excellently (within 7 %) with MLS measurement.R010 produces about 20-30 % less SO 2 near 26 km compared to R001, and significantly more above 30 km.This could be explained by the fact that R010 disperses slightly more SO 2 above 24 km compared to R001.The SO 2 vertical profile of R017 is shifted to lower altitudes compared with the observed values, likely due to its concentrated injection distribution near 19-20 km (see Fig. 1).Box14Mt and R153 fail to match the observed profile.SPARC20Mt agrees with the observations under 28 km better than Box14Mt and R153, but largely overestimates the observations above.The common feature of R153, Box14Mt and SPARC20Mt is that their initial vertical distributions release much more SO 2 above 24 km compared to R001, which is skewed towards lower altitudes, therefore retaining more than 90 % of emitted SO 2 below 24 km (Fig. 1).SO 2 profiles simulated by the two 3-D simulations (dashed curves in Fig. 2) are similar to the corresponding AER 2-D results, though SOCOL-AER predicts a lower maximum value and more readily distributes SO 2 to higher altitudes, reflecting differences in OH and transport between the two models.(Baran and Foot, 1994).SAGE II aerosol data is derived from the retrieval algorithm SAGE 4λ by Arfeuille et al. (2013), and includes only stratospheric aerosols.

Matching the burden
Figure 3 shows the evolution of the simulated stratospheric aerosol burden (megaton of H 2 SO 4 /H 2 O) compared to that derived from HIRS (Baran and Foot, 1994) and SAGE-4λ (Arfeuille et al., 2013).R001 matches the HIRSderived maximum aerosol burden of 21 Mt (equivalently 15-16 Mt of sulfate mass without water) during the first few months after the eruption, and after month 14 agrees with the SAGE-derived burden (mostly within 20 %).In contrast, SPARC20Mt reaches a maximum burden of 32 Mt of H 2 SO 4 /H 2 O, which is ∼ 50 % more than the 21 Mt derived from HIRS.R033 emits 17 Mt of SO 2 using the same vertical profile as R001, and peaks at 25 Mt of aerosol mass, about ∼ 30 % more than HIRS, whereas the uncertainty of HIRS is about 10 % (Baran and Foot, 1994).This means that the initial mass loading of 17 or 20 Mt of SO 2 into the stratosphere is apparently too high.Scenarios using 14 Mt of SO 2 show that the evolution of the aerosol burden is highly sensitive to different injection profiles.R010 initially distributes somewhat more SO 2 above 24 km compared to R001, and shows a better decay rate of the aerosol burden.R017 emits SO 2 mainly concentrated between 19-21 km, and its aerosol burden peaks similarly to R001, but declines more rapidly.R153 and Box14Mt inject about 60 and 40 % of their sulfur mass above 24 km, respectively, leading to a greater maximum aerosol burden and a slower decay rate of the burden than R001.R153 has even a slightly larger maximum aerosol burden than R033, though R033 has the larger initial SO 2 mass loading.Together, these results reveal that the injection altitude and initial mass loading affect the lifetime of the   (Deshler et al., 2003;Deshler, 2008), and model simulations in October 1991 (upper panels) and December 1991 (lower panels) for particle size channels r > 0.15 µm (left panels) and r > 0.5 µm (right panels).volcanic aerosol.An increase in the distance of the volcanic plume above the tropopause will increase the lifetime of the volcanic aerosol due to a longer residence time for sedimenting particles and a slower pathway of the aerosol within the Brewer-Dobson circulation.On the contrary, a larger initial mass loading may offset a higher injection altitude because of faster sedimentation caused by larger particles.
The results of "R001 3-D" using the coupled aerosolchemistry-climate model SOCOL-AER is consistent (mostly within 10 %) with the AER 2-D simulation R001.In contrast, the consistency between R153 and "R153 3-D" is less satisfactory.The maximum aerosol burden simulated by "R153 3-D" is within 10 % of R153, but the e-folding time of the aerosol burden in the 3-D simulation ("R153 3-D") is significantly faster (13 vs. 15 months) than in the 2-D simulation (R153).This indicates that in addition to the initial mass loading and microphysics, model dynamics is essential to the decay of the volcanic aerosols.This difference between R153 (AER) and "R153 3-D" (SOCOL-AER) is possibly due to an insufficient rate of exchange of air between the troposphere and stratosphere in the AER 2-D model (Weisenstein et al., 1997) and/or a faster Brewer-Dobson circulation with respect to observations in the SOCOL (see the "tape recorder" in Fig. 8 of Stenke et al., 2013).

Matching particle size distributions
Figure 4 shows comparisons between the optical particle counter (OPC) measurements operated above Laramie (Deshler et al., 2003;Deshler, 2008)  stations used within SAGE-4λ under conditions when the atmosphere was so opaque that SAGE II could not measure it; so-called gap-filled data with large uncertainty (SPARC, 2006;Luo, 2015).
0.5 µm).Below 23 km, R001 reasonably matches the observations for r > 0.15 µm, but less satisfactorily for r > 0.5 µm.The number density from R010 is slightly higher than R001 above ∼ 24 km, which is consistent with the comparison between initial vertical profiles of R001 and R010 (see Fig. 1).R017 agrees best with the observed number density, particularly above 24 km, because R017 emits very little SO 2 above 22 km.R033 predicts slightly higher number concentrations than R001 due to its larger initial mass loading (17 Mt SO 2 ), but shows in general similar results to R001.In contrast, the calculations from R153, Box14Mt and SPARC20Mt differ significantly from R001.Above 23 km, these three scenarios further overestimate the observations than R001 because their initial injection profiles release much more SO 2 above 23 km compared to R001.Be-low 23 km, R153 substantially underestimates the observations in October 1991 as its injected mass locates mainly between 23-27 km, while Box14Mt shows better agreement with the observations (r > 0.5 µm) below 18 km than R001, but largely underestimates the maximum near 21 km.SPARC20Mt is similar to R001 below 20 km since its initial mass loading (20 Mt SO 2 ) compensates for the deficiency of its vertical mass injection profile in the lower stratosphere.
The calculations from SOCOL-AER are generally consistent with the corresponding 2-D ones (R001 and R153).SOCOL-AER produces higher number concentration in October 1991 compared to the AER 2-D model.In December 1991, this difference between the 2-D and 3-D simulations shrinks, and "R001 3-D" further improves the agreement with the OPC measurements below 18 km for r > 0.5 µm.its higher initial mass loading, although it has the same vertical profile as R001.SPARC20Mt has even larger values than R033 due to a 20 Mt of SO 2 mass loading.Box14Mt and R153 largely overestimate the observed extinctions above 24 km.The 3-D simulation "R001 3-D" is superior to all the 2-D simulations, while "R153 3-D" performs worse than the 2-D simulations R001 and R033.Likewise, in June 1992, R001, R010 and R017 also do a better job than other 2-D simulations.The two 3-D simulations "R001 3-D" and "R153 3-D" are now both superior to all 2-D model results, although the differences between them start to shrink as the their aerosol burdens are now within 10 % from each other.
Here the 3-D model shows a better extinction vertical profile likely because the 3-D model uses an improved numerical scheme based on Walcek (2000) for sedimentation, while the 2-D model uses an upwind scheme, which would cause artificial upward transport of particles to high altitudes (Benduhn and Lawrence, 2013;Sheng et al., 2015).Evaporation of aerosol becomes important only above ∼ 32 km, therefore should play a minor role in explaining the discrepancy between the 2-D AER and 3-D SOCOL-AER.Overall, the results from SPARC20Mt, Box14Mt, R033 and R153 display a common deficiency, as they tend to overestimate aerosol extinctions in high altitudes above 24 km.Excessive mass loading (as in SPARC20Mt or R033) is one of the reasons.However, the shape of the initial mass vertical profiles appears to be at least as important as the initial mass loading.Box14Mt has 30 % less total mass loading than SPARC20Mt, but it shows even higher extinctions in high altitudes because it has 40 % of its mass injected above 24 km, while SPARC20Mt has only about 20 % of its mass there.
Figure 6 compares the modeled aerosol optical thickness (AOT) with the SAGE II measurements.The southward transport of the volcanic clouds observed by SAGE II is reasonably reproduced by the models.The best scenarios here are R001 and R010, whose SO 2 injection profiles peak between 18-21 km and disperse the volcanic plume broadly (σ = 4 km).In contrast, R017 with a narrow dispersion (σ = 2 km) constricts the initial SO 2 between 18-22 km, which leads to a faster decay of AOT than R001 and R010.R153 and SPARC20Mt distribute too many aerosols to high latitudes due to injecting SO 2 excessively above 24 km.The impact of the initial SO 2 vertical profile on the hemispheric dispersion of the volcanic clouds is more pronounced in the 3-D simulations as shown in the two bottom panels.These results show that spatiotemporal distribution of the volcanic aerosols is affected by initial injection profile of SO 2 and the optimal parameters found in Table 1 would lead to better model results when compared to SAGE II observations.

Conclusions
We have conducted over 300 Pinatubo-like simulations by perturbing four parameters which determine the magnitude and vertical distribution of injected SO 2 .Our simulations show that the initial SO 2 magnitude and distribution play a significant role in the evolution of stratospheric aerosol properties following the Pinatubo eruption, including rates of poleward transport of volcanic clouds.Our ensemble study suggests that Pinatubo injected less than 17 Mt of SO 2 into the stratosphere, and that good agreement can be reached with a 14 Mt injection.The vertical profile of the injected SO 2 is likely skewed towards the lower stratosphere, with 80 % of the SO 2 mass injected below 24 km and the distribution peak likely between 18 and 21 km.We have found a set of initial injection parameters such that the resulting model simulations fairly reproduce the evolution of stratospheric aerosol properties when compared to HIRS and SAGE II based data.This reduces the uncertainties in modeling the initial sulfur mass loading of Pinatubo.

Figure 1 .
Figure 1.Vertical distribution function F (z). Black line: used in SPARC (2006).Blue line: uniform (box) profile that distributes SO 2 homogeneously with altitudes.Each of these curves encloses a unit area.

Figure 3 .
Figure 3. Evolution of simulated global stratospheric aerosol burden (Mt H 2 SO 4 /H 2 O) compared to the HIRS and SAGE II-derived data.HIRS-derived data include both tropospheric and stratospheric aerosols(Baran and Foot, 1994).SAGE II aerosol data is derived from the retrieval algorithm SAGE 4λ byArfeuille et al. (2013), and includes only stratospheric aerosols.

Figure 5 .
Figure4shows comparisons between the optical particle counter (OPC) measurements operated above Laramie(Deshler et al., 2003;Deshler, 2008) and model-calculated cumulative particle number concentrations in October and December 1991 for two size channels (r > 0.15 µm and r >

Table 1 .
Scores and rankings of 326 AER 2-D atmospheric simulations of the Pinatubo eruption sorted according to the weighted rank ("RankWt").The weighting is given by 16.7 % of the SO 2 score (ScoreSO 2 ), 16.7 % of the OPC score (ScoreOPC), 33.3 % of the global burden score (ScoreBurden), and 33.3 % of the aerosol extinction score (ScoreExt).The rank computed by the arithmetic average of the four scores is also provided ("RankAvg").Scores of two additional 3-D simulations "R001 3-D" and "R153 3-D" from the aerosol-chemistryclimate model SOCOL-AER are provided at the bottom of the table.