Radiative forcing by volcanic eruptions since 1990, calculated with a chemistry-climate model and a new emission inventory based on vertically resolved satellite measurements

This paper focuses on the injection of SO2 to the stratosphere by volcanic eruptions, and the resulting variability of the stratospheric aerosol layer. It presents a new volcanic SO2 emission database, derived from a collection of satellite instruments, covering the period 1990-2019. It also presents results from a chemistry climate model which uses the updated injection database, and compares the results of the model to various satellite data sets, focusing on the multi-wavelength aerosol optical depth and instantaneous radiative forcing produced by the aerosols.

The construction of such detailed SO2 injection estimates covering the 1990-2019 period is an impressive accomplishment. It is also to my knowledge quite novel, as I believe it is the first attempt to produce SO2 injection values from sulfate aerosol extinction measurements. Unfortunately, the description of the methods used to produce these estimates is lacking. Furthermore, assumptions and choices made in the methods are not given justification. More detailed comments are included in "Major Comments" below.
Chemistry climate model simulations using the new SO2 injection data set are performed and some results shown. Good agreement with observations is achieved, but there is insufficient analysis to provide any improved understanding of the physical or chemical processes that control stratospheric aerosol evolution.

Major comments:
The description of how SO2 amounts were calculated lacks sufficient detail. I am not aware of any other study that has estimated SO2 injection amounts based on aerosol extinction measurements. This is thus a novel technique, but the method used is not described beyond a few statements along the lines of "The SO2 mixing ratio perturbation is derived from the extinction perturbation observed in a 10-day period beginning about a week after the eruption by dividing by air density, multiplying by a constant and subtracting a typical background." This explains extremely little: what constant is used, and why? How is the typical background determined? How well can the volume of the aerosol cloud be estimated a week after eruption from the satellite measurements? SAGE in particular has a very sparse sampling density, how does this impact the estimates? Can the method be validated? It would seem that the method could be applied to SAGE and OSIRIS during periods of overlap with MIPAS and the values from the new method compared to the "direct" MIPAS measurements. This would help increase confidence in the method, and provide some idea of the uncertainties in the estimates. -> Response: We added a detailed explanation with case studies in Appendix B: "The eruption of Reventador in the tropics in November 2002 has shown to be an ideal case where simultaneous observations of all satellite sensors were available so that the direct SO2 observation could be used for development and validation of a conversion formula for the 750 nm extinction seen by GOMOS and OSIRIS, which works also approximately for SAGE if its observations at 530 and 1025 nm are interpolated to 750 nm. Here we first use the ratio between model calculated sulfate volume mixing ratio and its share on extinction in low latitudes of the lower stratosphere which is typically 1.2×10 12 ×air density (in molecules/cm 3 ). This works for medium size eruptions and data available over about four weeks following the eruptions, and if no other events occur less than about four weeks before which is the case for the Reventador eruption. If the time lag of data is several weeks a correction factor >1 has to be applied to account for removal processes, if another event is relatively close in time, the factor has to be <1 to remove the influence of the previous event. For Reventador the factor is 1 (for OSIRIS 0.8 is slightly better). From all instruments the derived injected SO2 mass is very close to 77 kt as shown in Table 2. The spatial patterns are similar, except when the zonal wind causes a shift in longitude due to the time lag from conversion of SO2 to aerosol, see Figure B1. In the case of SAGE, the alternate method of Grainger et al. (1995) involving aerosol surface area density (SAD) and aerosol volume density is more suitable to remove cloud perturbation. It is assumed that sulfate mixing ratios correspond to the SO2 injected. Some uncertainty remains from removing the background which we have done by subtracting a fraction of the derived SO2 at the longitude where it has a minimum, i.e. the longitude where the effect of the volcano is smallest for all altitudes. Integrated injected SO2 masses for all examples are provided in Table B1. For the eruption of Merapi in November 2010 the satellite instruments do not agree. From OSIRIS about 70% more injected SO2 is derived than from MIPAS, i.e. 170 kt instead of 97 kt used in the transient simulation (see Table 2 and differences in Figure 10). GOMOS has too sparse data here to obtain a proper integral directly but patterns are similar ( Figure B2). If other information is available, the gaps can be filled with likely values in the region where the plume was seen, a method which had to be applied also to some events seen by OSIRIS in 2018 and 2019 for which the data were sparse. For high latitude eruptions the longer conversion time of SO2 to sulfate compared to the tropics has to be considered which, together with aerosol removal processes, lead to a weaker extinction signal. To account for this a correction factor of about two in the conversion formula for OSIRIS for example for Sarychev in June 2009 leads to values consistent to the ones derived by MIPAS ( Figure B3). For the low latitude eruption of Mando Hararo in the same entry of Table 2 (separated at 24° N for the integration) the factor 1 is still appropriate." I highly recommend that the emission database be provided as an electronic supplement (e.g., csv or xls), to allow it to be readily used by other researchers. -> Response: The input data files and model output of EMAC used here are stored at DKRZ, Hamburg, the volcanic inventory and the output for radiative forcing also at WDCC https://doi.org/10.26050/WDCC/SSIRC_3 The table, as text, presently takes up almost 8 pages of the manuscript: it would be more efficient to visualize the data somehow and include the values as supplemental information. Also, I strongly suggest that the format of the table be modified so that each individual eruption be listed per row, even if there are multiple eruptions on a given date. This will greatly improve the ease in which the data can be read within a computer program and thus used in other studies. -> Response: SO2 mixing ratios from the volcanic emission inventory are shown in Fig. 6. As an essential part of the novelty of this paper, the table should remain in the text because it is a comprehensive reference that cannot be represented by a single visualization. Additionally, the table is available for data processing as formatted ascii at this link: https://doi.org/10.26050/WDCC/SSIRC_3 The model results show good agreement with observations, but it's impossible to know whether the improved agreement (compared to prior works from the same group) is a result of the updated SO2 injection data, or to model improvements or changes in model resolution. Given the theme of the ACP journal, the reader expects that this work should improve our understanding of the chemical and/or physical processes that control stratospheric aerosol evolution, but it remains unclear if there is any improvement in understanding being extracted from the study. Nor is there any real motivation or objectives stated in the introduction for the model simulations.
-> Response: The improvements are mostly due to use of more satellite data for volcanic SO2 using a novel method. This is mentioned now at several places including the abstract. Model improvements include consideration of aerosol effects in the photolysis rate calculation (minor effects here) and compared to Bingen et al. (2017) a finer horizontal resolution (see section 3).

Specific comments:
L11: "Reproduce" is too strong -> Corrected: are consistent with L12: Here it is said that "slight deviations … were found only for the large volcanic eruption of Pinatubo in 1991", but later in the document deviations in other time periods, e.g., 2010 are discussed, so this is inconsistent. -> Removed: only L19: precise language is needed here, is this the peak radiative forcing produced by a typical "small" eruption, or the time average forcing from these eruptions? And what is a small or medium eruption? Also, it's not clear how this number is estimated, a value of 0.10 W/m^2 is not mentioned in the results or conclusions, and if it comes from Fig 11, how is the effect of small eruptions separated from that of "background" sulfur (e.g., DMS, OCS) transported into the stratosphere via atmospheric circulation? -> Response: This number can be taken from Fig. 11 for volcanically quiescent periods or periods between medium size eruptions. Background is about 0.04 W/m 2 (not shown explicitly, taken from a simulation with much less volcanoes). This includes organics and dust.
L22-24: references needed for these statements. L37: I am skeptical of a 3-year upper limit on the impact from volcanic eruptions: if ocean temperatures are a part of "climate", then there is good evidence that volcanic impacts on climate can last much longer than 3 years (e.g., McGregor et al., 2015). Obviously the period of impact depends on many factors, but we should be careful to not overly simplify statements which might be misleading to some readers. L65: Some information should be given on how the SO2 column data was used, especially in regards to how a stratospheric component was estimated from the full column. -> Response: In the case of data gaps for the four main satellite instruments used, the data from additional satellites e.g. TOMS, OMI or OMPS are applied to double check, if available. Especially in 2018 and 2019 the OSIRIS data are so sparse that constraints from instruments like OMPS or analogues events of previous years have to be superimposed for some eruptions.
L130: The gaps in spatial coverage of the OSIRIS data at 17 km extend significantly beyond the polar night: they seem to extend even in best cases to 20-30deg. Some rephrasing needed.
-> Response: OSIRIS provides a surface coverage from 82° S--82° N, except in polar winter when there is no sunlight and except in the Southern Hemisphere winter when tangent point is not illuminated by the sun L136: It's not apparent how the sensitivity to clouds can be seen in Figure 4. -> Response: At altitudes near and below the tropopause, the OSIRIS measurements are sensitive to clouds that may be interpreted as elevated aerosols. This is likely contributing to larger background extinction values measured below approximately 17 km in the tropics, as can be seen in Figure 4. L140: How is the correction factor determined? This sounds suspiciously like numbers have been chosen only to produce best agreement. -> Response: To estimate the factor, we iterated calculated extinctions to agree with OSIRIS and also used observations and assumptions by Vernier et al. (2016). A detailed description with 3 case studies is added to the Appendix: "If the time lag of data is several weeks a correction factor >1 has to be applied to account for removal processes, if another event is relatively close in time, the factor has to be <1 to remove the influence of the previous event." L157: The study of Grainger et al. (1995) does not seem to provide a relationship between SAD and SO2 mixing ratio. More explanation needed. -> Response: This is skipped here. Now we write in Section 4: "For SAGE II in most cases the SO2 mixing ratio is derived using the parameterisation of Grainger et al. (1995) which converts SAD to volume density as a first step. We use the pressure and temperature provided to convert from mass density to a volume mixing ratio, assuming that observed sulfate is produced from injected SO2 some weeks ago. With this method it is easier to correct for cloud contamination than by using the extinction directly as above for the other instruments. Case studies for three events, comparing SO2 results from the different satellites and the different conversion methods are presented in Appendix B".
L190: It is not clear how differences in the "vertical transport of tracers, like dust and water vapor or ozone" between model resolutions has any importance to the present study.
-> Response: This is part of the general setting of the model simulations and has been moved to the appendix.
L216: What parameters? -> Response: We removed "parameters" by: "aerosol optical properties like wavelength-dependent particle extinction cross section, single scattering albedo, and asymmetry parameter for each aerosol mode from AEROPT (Dietmüller et al. 2016) " L218ff: The double radiation call most likely calculates the "instantaneous radiative forcing". It is important to be clear about this and consistent with the terminology. -> Added: "instantaneous" radiative forcing L219: There is a double radiation call, but how exactly is the radiative forcing calculated? -> Correction: "instantaneous" radiative forcing, "..by taking the difference of the net total fluxes at 100 hPa or TOA." L220: Not understanding this, are you diagnosing the impact of volcanic aerosol on upper stratospheric UV absorption? Nothing like this is shown in the results. -> Response: The description of the RAD\FUBRAD sub-submodel is part of the general setting of the model simulations and has been moved to the appendix. L241: What is the justification for the lower limits to the vertical integration given? You use 12 km as the lower limit in high latitudes, but the climatological tropopause height in high latitudes is 9-10 km. Conversely, you use 14 km in low latitudes, but the tropopause there is around 17 km. A thorough explanation for these counterintuitive thresholds will need to be given. -> Response: The lower limit of 12 km altitude at high latitudes was chosen based on the signal-tonoise ratio, uncertainties for low altitudes, and clouds in the volume mixing ratio profiles obtained by MIPAS and the other used satellite instruments (e.g. effects of frontal systems). In the tropics, we set the lower limit at 14 km to account for transport processes in the UTLS layer, especially during the Asian summer monsoon. Here in our extraction scheme we exclude cloud contaminated regions (see section 4 and 2). L251: An "integration time" has not been introduced, it is not clear what this means in terms of the method.
-> Response: The temporal resolution of the satellite data is 5 days (for MIPAS, GOMOS and OSIRIS). The integration time is the case dependent time period used for the single eruptions. "The integration time (i.e the used time period)" Table 2: There are a number of cases where the number of values do not match between the different columns in a particular row, e.g., 11 Feb 1990, 19 Aug 1992, 18 Sep 1996. Expanding the table so each eruption is listed in a single row would help this issue, as well as improve the machine readability of the table more generally. There is also a case (14 Jan 2002) where values are listed within brackets, and I did not find an explanation for what this means. -> Response: These issues have been corrected, the table is available now in a Fortran formatted form and it has to be consistent to the 3D netcdf SO2 pertubation files provided in https://doi.org/10.26050/WDCC/SSIRC_3, which often contains multiple events. A single event list with more than 500 entries would not be acceptable for the ACP readership.  Guo et al., 2004), but in contrast to recent model studies which suggest the effective injection for Pinatubo was much less (e.g., Mills et al., 2016;Dhomse et al., 2014). Some discussion of this issue would fit well into the paper. -> Added in sec. 6.2: "In this study about 17 Tg SO2 are injected for the Pinatubo eruption (Guo et al. 2004 2013)). Thus, this study is in the middle range of the injected sulfur mass. On the other hand, filling the gaps in the SAGE data just by horizontal linear interpolation increases the peak AOD by about a factor of 2, which is close to the GloSSAC compilation. In Figure 10 the AVHRR (Advanced Very High Resolution Radiometer) by Long and Stowe (1994) at 630 nm are included, which are close to our simulations (if converted to 750 nm their AOD would be slightly less or to 550 nm slightly larger). When comparing the EMAC simulations (red line in Figure 9) with the simulation of Schmidt et al. (2018, fig.1) (black line in Figure 9, lower panel) it can be recognized that a smaller value for the peak of the Pinatubo eruption occurs, but here it needs to be considered that Schmidt et al. (2018) are using monthly global-means. This has the consequence that the signal of single eruptions is blurred and smaller sized eruptions cannot be easily identified."  L293: Is the statement on SO2 lifetimes made here a result of this study, or are the lifetimes equivalent to those given by Hoepfner et al. (2015)? If the result is the same as Hoepfner et al., (2015), that should be explicitly stated. If estimated lifetime are different from Hopefner et al. (2015), how and why? -> Response: "Generally, the conversion of SO2 to sulfate aerosol particles depends on several factors, such as the altitude, latitude, or season of the eruption and takes according to Höpfner et al. (2015) about 13, 23 and 32 days in 10-14, 14-18 and 18-22 km altitude, respectively, in midlatitudes. Carn et al. (2016) report an e-folding time varying between 2-40 days. The range agrees with our simulations (and assumptions in section 4)." L300: This sentence seems to say that stratospheric aerosol optical properties were calculated using a range of different aerosol types (sulfate, dust etc.). Is this correct, or is the sentence just misleading? -> Response: Yes, this is correct. See Dietmüller et al. (2016): "Aerosol species explicitly considered are water soluble inorganic ions (WASO), black carbon (BC), organic carbon (OC), sea salt (SS), mineral dust (DU), and aerosol water (H2O). The refractive indices for those aerosol species are extracted from various data sources (most of the data are compiled in the HITRAN2004 database) and include wavelength dependencies. … The refractive indices for each aerosol mode required as input for the lookup tables are calculated assuming an internal mixture of the aerosol components for the hydrophilic modes. A mean refractive index is calculated for each mode wavelength combination by averaging the refractive indices of the individual components weighted with their volume contributions. The corresponding Mie size parameters are derived from the median radii of the log-normally distributed modes and the respective wavelengths. The wavelength-dependent particle extinction cross section, single scattering albedo, and asymmetry parameter for each mode are then obtained from the lookup table for the appropriate modal width (σ)." L330: The OSIRIS data is converted from 750 nm to 550 nm, which is fine, but this contradicts the statement just a couple sentences earlier that "Unlike most other studies, the stratospheric AOD is compared at the original wavelengths derived from different optical channels of the satellite instrument measurements." -> Response: The model calculates the AOD and other optical properties directly, derived from the original wavelengths of the satellite data. For additional comparisons only, the satellite data of OSIRIS were converted to 550 nm by the cited authors.
L333: The statement that "differences after the large Pinatubo eruption in 1991 between the model simulations and the SAGE II observations are related to the "saturation" effects of the satellite instrument" seems much too confidently worded. It seems quite possible that "saturation" effects explain some of the difference, but how certain can you be sure that it is the only, or even the primary reason? In the tropics, the simulated AOD appears to be ~3 times larger than the SAGE II measurements-is it likely that the SAGE II measurement is so strong an underestimate of the true total AOD? -> Response: See AVHRR data now included. Filling the SAGE data gaps in the lower stratosphere just by linear interpolation in the first year after the eruption increases AOD by about a factor of 2. L362: clarify that the *simulated* AOD drops too quickly compared to the observations. -> Corrected: "the simulated instantaneous global negative radiative forcing drops again too quickly for the EMAC simulations (red line) after the Pinatubo eruption compared to the observations." L386: This appears to be a result of the study by Bruehl et al. (2018), which would be important in describing the experiment earlier in the manuscript but not here in the conclusions. -> Response: Skipped here (see above).
L388: The SAGE II and OSIRIS extinction measurements are not really "newly available", some version of this data has been available for many years. The estimation of SO2 from these data sets is quite new-it's what this paper is presenting! -> Response: The resolution of the updated versions is improved, so OSIRIS data has allowed comparisons up to be extended to 2019; additional SAGE II data was also used extend the comparison back to 1990 (together with the Smithsonian database).
L402ff: This conclusion is not supported by the results: there is no quantification of the impact the increased number of eruptions included in the database has on the radiative forcing, or its level of agreement with observations. -> Removed.
L408ff: This is an interesting conclusion, but it is not supported by the results. There is no demonstration that including the injections below the tropical tropopause improves the agreement. Even a comparison with prior studies will not prove necessarily support the statement since those prior studies used a different resolution model. -> Response: Remarks on that are added at several places in other sections.
L418: This is not a new result, as it has been shown by prior studies. -> Response: Paragraph skipped. Instead we included now in the paragraph beginning with: "Our volcanic sulfur emission inventory…." "The inclusion of plenty of small size eruptions reaching the UTLS has the consequence that stratospheric aerosol optical depth and radiative forcing does not decrease to almost zero between medium size eruptions in agreement with observations, in contrast to a lot of other studies." L422: The impact of volcanic aerosol on tropical upwelling is not diagnosed in this study. Prior studies have explored this, but statements like this can not be included in the conclusions of this work if there are no new results shown to support it and build upon prior work.
L437ff: This paragraph talks about meteoritic dust, which was not investigated in the study. Perhaps simply adding a sentence or two on the agreement between the model and observed aerosol extinction in the upper stratosphere to motivate the discussion of meteoritic dust would help the reader follow the logic here. -> Removed: Remark added to sec.6.1: "Above about 24 km altitude, EMAC underestimates the observations because in the model meteoric dust particles were not considered." 448: Confirming the findings of the IPCC report is, firstly, incorrectly phrased, since the IPCC report only summarizes and reports findings gathered from the published literature. It would be more important to compare the results here with the primary sources, including studies that have been published since the IPCC AR5 (e.g., Schmidt et al., 2018). Second, confirming some general results from prior studies does not make an overwhelming case for publication. What does this study add to the understanding of volcanic radiative forcing that wasn't known before? -> Response: Paragraph shortened, restricted on new results. L450: Radiative forcing is stated to be that at the surface here, where Fig 11 is said to be RF at the tropopause. Also the numbers quoted here don't seem to agree with Fig 11. It would be best to only refer to calculations for which the results are shown in the paper. -> Response: Typo corrected, numbers refer to difference to quiescent periods which is now mentioned.

Editorial comments:
Line 9: Volcanic SO2 is not "pollution" in the usual sense of the word, suggest it be cut here. -> Corrected.
L55: Awkwardly phrased: the processes aren't structured, the paper is, and not strictly according to processes. -> Corrected: "This paper is structured as follows:" L80: I've never seen pptv written with v as a subscript, is this a new standard? -> Corrected.