Interactive comment on “ A perturbed parameter model ensemble to investigate 1991 Mt Pinatubo ’ s initial sulfur mass emission ”

The conclusion that the lowermost estimate of earlier studies for SO2 injected by Pinatubo should be selected is only of limited value because it appears to be significantly perturbed by model artifacts and arbitrary weighting of ’scores’ for differences to observations. Giving the MLS SO2 measurements a higher weight and the rather uncertain burdens estimated from SAGE during saturation of the instrument a lower one would completely change the conclusions. From text and figures it is also often not clear how the scores were calculated. In the table and the figures important cases are missing.


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 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.

30
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).

35
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 40 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 measure-45 ments 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 three months after the eruption.Given the lack of measurements in the pe-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 quasibiennial 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 Weisenstein et al. (1996) :::::::::::::::: Notholt et al. (2005).Photodissoci-110 ation 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 repre-115 senting 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 humid-120 ity 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, thus 125 the entire particle size distribution.Operator splitting methods are used in the model with a time step of one hour for transport, chemistry, and microphysics, and 3-minute substeps for the microphysical processes that exchange gasphase H 2 SO 4 with condensed phase, and 15-minute substeps 130 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)).

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 into 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 (e-folding time) is sufficiently slow (τ ∼ 25 days) and the zonal transport around the globe sufficiently fast (τ ∼ 20 days) (Guo et al., 2004) 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 the satellite extinction 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 Limb Sounder (MLS) measurements (Read et al., 1993).
observations.Similarly, skewness α > 0 can lead 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 A selected list from the 326 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 225 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 R149 3-D at the bottom of Table 1) using the coupled aerosol-chemistry-climate model SOCOL-AER 230 (Sheng et al., 2014) :::::::::::::::: 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 R149) due to different vertical model levels between the two models.
Metrics and data sets.

320
In contrast, SAGE II, as an occultation instrument, becomes very reliable when the stratosphere starts to be sufficiently transparent.The measurement uncertainty is generally better than ∼20% for 525 nm wavelength and ∼10% for 1020 nm (see Fig. 4.1 in SPARC ( 2006)).Therefore, ScoreExt is 325 weighted as one third for 525 nm and two thirds for 1020 nm.Finally, ScoreBurden uses the HIRS-derived data up to month 12 :::::::: 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 330 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 335 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 re-  (Baran and Foot, 1994).SAGE II aerosol data is derived from the retrieval algorithm SAGE 4λ by Arfeuille et al. ( 2013), and include only stratospheric aerosols.
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 year after Pinatubo (see Figure 3).After this period, HIRS tends to be noisy due to its lack of sensitivity at high latitudes where there is a contribution from errors in the background signal (Baran and Foot, 1994).In contrast, SAGE II, as an occultation instrument, becomes more reliable when the stratosphere starts to be sufficiently transparent.Therefore, ScoreBurden uses the HIRSderived data up to month 12 :::::::: December ::::: 1991 and the SAGEderived data afterwards, with an overall uncertainty of 20%.
Table 1 shows the scores of selected scenarios, sorted according to the weighted rank ("RankWt" in the next to last column).The best scenarios (RankWt ≤15) 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 19-22 ::::: 18-21 km with 3-4 km width (scale parameter σ).Location parameters µ larger than 22 km are generally 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 Figure 1.The worst scenarios (RankWt ≥317) 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 R149 rank much more poorly than the optimal scenarios, although their injection mass is the same, because their vertical profiles (shown in Figure 1) inject over 50% mass above 23-24 km.The sce-375 nario R034 has the same vertical profile as R001, but more emitted mass (17 Mt SO 2 ), leading to poorer ranks in the aerosol burden and extinctions.The scenario SPARC20Mt ranks at 211 in Table 1, although its vertical profile is close to the optimal scenarios (about 10-20% more mass above 23 380 km).This implies that emitting 17 or 20 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) 385 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) 390 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).
The common feature of R149, 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 (Figure 1).SO 2 distributions :::::: profiles : simulated by the two 3-D simulations (dashed curves in Figure 2) are similar to the corresponding AER 2-D simulations ::::: 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.Matching the burden.
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.Below 23 km, R149 substantially underestimates the observations in October 1991 as its injected mass is mainly situated :::::: 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 calculations :::: ones (R001 and R149).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.
We compare the modeled 1020 nm extinctions with the gapfilled SAGE II version 7.0 (Figure 5).SAGE II data points with horizontal bars are actual SAGE II measurements and denote natural variabilities, while data points without bars are gap-filled from lidar ground stations, which have a higher uncertainty (SPARC, 2006).Figure 5 shows comparisons in January (upper panel) and July (lower pannel) 1992 for five latitude bands from left to right: 50-40 • S, 30-20 • S, 5 • S-5 • N, 20-30 • N and 40-50 • N.

Conclusions
We have conducted over 300 Pinatubo-like simulations based on variations of four parameters of initial total SO 2 mass and altitude distribution.These parameters control the temporal and spatial evolution of stratospheric aerosols in the years following the Pinatubo eruption.The altitude distribution of SO 2 injection is represented by a skew-normal distribution.
Our simulations suggest 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, 80% of which was injected below 24 km with the maximum located between 19-22 ::::: likely ::::::: between ::::: 18-21 : km.This reproduces HIRS and SAGE II-based estimates of the evolution of total stratospheric aerosol burden.Furthermore, this largely improves the previous overestimates in :::::::: presented :: in :::::::::::::: SPARC ( 2006) in : modeled extinctions at high altitudes when comparing to SAGE II gap-filled measurementsSPARC ( 2006), and realistically simulates aerosol extinctions in the lower stratosphere.We have defined an optimal set of the emission parameters such that the resulting burdens and extinctions match satellite and lidar measurements, and reduce the uncertainties in modeling the initial sulfur mass loading of Pinabuto.
') 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 Figure 1); (2) SPARC20Mt is the reproduction of the Pinatubo simulation conducted in SPARC 220 (2006), which injects 20 Mt of SO 2 and has a vertical profile 'SPARC' shown in Figure 1.

Fig. 3 .
Fig.3.Evolution of simulated global stratospheric aerosol burden (Mt H2SO4/H2O) 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 include only stratospheric aerosols.

Table 1 .
Scores and rankings of 326 AER 2-D atmospheric simulations of the Pinatubo eruption sorted according to the weighted rank ("RankWt"). ::: table.