Interactive comment on “ Retrieval of volcanic SO 2 from HIRS / 2 using optimal estimation ”

The manuscript presents a new scheme for the retrieval of atmospheric SO2 total column amounts after volcanic eruptions from HIRS/2 observations. An optimal estimation retrieval using radiances from three HIRS/2 channels in the mid-infrared region is presented. Retrieved parameters are cloud top and the total column amounts of water vapour and SO2. Major error sources, which are identified by synthetic observations, are cloud/ash interference and the assumptions on the altitude and vertical extend of the SO2 plume. Through simulations it is further demonstrated that the new scheme is superior compared to simple brightness difference methods. The method has been applied to the case of the eruption of the Cerro Hudson volcano in 1991.


Introduction
Volcanic eruptions are important for climate and climate change.They perturb atmospheric chemistry and radiative transfer.Their signal in climatic records must be accurately quantified before any attribution of climate change to anthropogenic sources.Furthermore, by studying the response of the atmosphere to volcanic eruptions in terms of climate sensitivity, one can test ideas relating to climate prediction.
The monitoring of volcanic SO 2 emissions, the main precursor to sulfate aerosols, is crucial not only for accurately characterising total emission estimates but also for understanding plume evolution.Until the mid-1990s, only one principal instrument (the Total Ozone Mapping Spectrometer, TOMS) has been able to observe eruptions for an adequate period to generate something approaching a climate relevant record.The sensitivity of TOMS limits it to detecting only the larger, explosive eruptions rather than effusive ones where material remains predominantly in the troposphere.Satellite instruments that have been used to measure volcanic SO 2 are given in Table 1.From 1996, with the ad-Table 1. Instruments (many of which were flown aboard several different platforms which are not listed) that have been used to measure volcanic SO 2 in the atmosphere.

Instrument name
Viewing geometry, Period of Relevant reference spectral region operation TOMS, TOMS-like instruments (e.g.SBUV/2) Nadir, UV 1979+ Krueger (1983), Kerr et al. (1980); Krueger et al. (1995Krueger et al. ( , 2008)); Guo et al. (2004 vent of the Global Ozone Monitoring Experiment (GOME) class instruments (UV-vis spectrometers), sufficient spectral resolution (and spatial resolution) has enabled the detection of lower amounts of SO 2 with higher accuracy from increasingly smaller eruptions.This has improved further still with instruments such as the Infrared Atmospheric Sounding Interferometer (IASI), from which SO 2 , sulfate aerosol and ash may be derived simultaneously due to its high spectral resolution and broad spectral coverage (Karagulian et al., 2010).Total erupted mass estimates for volcanic eruptions can often differ by greater than 100 % between instruments, as a result of sampling; geometry; differences in sensitivity; and assumptions that contribute to algorithms, such as plume height.For example, Thomas et al. (2009) present a multi-sensor comparison of the 2005 eruption of Sierra Negra (Galapagos Islands), using concomitant observations by TOMS, OMI and MODIS.They found a wide estimate of total erupted SO 2 calculated from the three instruments, ranging from 60 to 1800 kT.It is still the case that the operational period of these more sensitive, recent instruments is not yet long enough to constitute a climate-relevant record.Here we present the methodology for a relatively fast and accurate volcanic SO 2 detection and quantification method for an instrument originally designed to operationally measure water vapour and temperature profiles.
The High-Resolution Infrared Radiation Sounder 2 (HIRS/2) has the potential to have captured stratospheric emissions from explosive eruptions continuously since 1979, and with significantly higher temporal sampling and greater sensitivity than TOMS.This enables the 35-year volcanic SO 2 emission record from satellites to be significantly enhanced, with potential uses for constraining models and examining in detail individual eruptions and plume evolution.
1.1 HIRS/2 instrument HIRS/2 is one of three instruments that originally constituted the Television Infrared Observation Satellite (TIrOS) Operational Vertical Sounder (TOVS), designed to provide atmospheric profile measurements of temperature and water vapour structure (Smith et al., 1979).The other TOVS instruments were the Stratospheric Sounding Unit (a radiometer) and the Microwave Sounding Unit (a scanning microwave spectrometer).The TOVS suite of instruments was first launched in 1979 aboard the new NOAA satellites based on the TIrOS-N (first of its class) design and evolved into the Advanced TOVS (ATOVS) system.Subsequent replacements have been deployed for the last 30 years aboard NOAA satellites (NOAA 6-17) (JPL, 2003) and more recently European platforms including most recently MetOp-A and MetOp-B as HIRS/4.Throughout its deployment there have been at least two instruments (and occasionally three) orbiting simultaneously.HIRS/2 has 19 detector channels in the infrared and one in the visible part of the spectrum for cloud detection during the day.These channels are relatively broad, spanning between 0.1 and 0.5 µm depending upon wavelength.The key instrument parameters are given in Table 2. Two HIRS/2 channels coincide with SO 2 spectral absorption features, these being 7.3 µm (a strong asymmetric stretch vibration band) and 8.6 µm.The precise central wave number is dependent upon instrument version, and only HIRS aboard NOAAs 10 and 12 featured an 8.6 µm channel.These channels were originally chosen to be sensitive to water vapour for use in sounding and applying corrections for the CO 2 and window channels.The 8.6 µm channel is also reported to be sensitive to volcanic ash and other aerosols (Kearney and Watson, 2009).
Channel 11 from HIRS/2 aboard NOAA 11, centred on 7.2 µm, is shown in Fig. 1.Also shown are simulated transmission spectra for water vapour (which this channel was designed to detect) and SO 2 , for two column amounts (1 and 300 DU).It demonstrates both that the channel and spectral feature coincide well and that for large column amounts of SO 2 the channel would be strongly affected.
1.2 Previous efforts to retrieve of SO 2 with the HIRS instrument Prata et al. (2003) demonstrated a method to detect volcanic SO 2 from HIRS, providing the SO 2 perturbation is strong enough and located above any significant sources of water vapour.It is based on a synthesis of the expected clean atmosphere brightness temperature (BT) for the channel and the observed deviation from it when contaminated by SO 2 .This method, hereafter referred to as either the Prata fit method or after Prata et al. (2003), uses a linear interpolation between the brightness temperatures of adjacent channels.It also assumes a fixed height of erupted volcanic SO 2 , since theoretically only one piece of information can be obtained from one channel, and column amount is not insensitive to the height of the plume.The technique requires the SO 2 to be located in the upper troposphere-stratosphere above most of the atmospheric water vapour, and there is no information about 17.4 km diameter Ground resolution IFOV (end of scan) 58.5 km by 29.9 km Distance between IFOV 42 km along track and nadir the height of the plume from the instrument itself.This information may be gleaned from other types of observations, but the fit is reliant upon the accuracy of this independent information.A description of how the Prata method operates is detailed in Prata et al. (2003).Whilst useful in itself, its most significant shortcoming is that, due to its simplicity, the model is unable to capture atmospheric variability (other than potentially that of SO 2 ).This particularly alludes to the variability of cloud, temperature and water vapour.Without independent height information of the SO 2 the radiance relationships are subject to potentially significant error.Indeed, it is not possible to formally quantify error of mass estimates from this method as it currently stands.Its strengths are that the operations required are computationally inexpensive and straightforward, as it is based on the principles of a band model.It has also performed well against other observational data sets, although the previously mentioned uncertainties that contribute to error make quantifying overall uncertainty difficult.It uses a minimum offset threshold in brightness temperature for the channel affected by SO 2 in order to predict the presence of SO 2 and yet excludes the effects of atmospheric water vapour variability.As such, its sensitivity to low amounts of SO 2 is limited.Guo et al. (2004) presented a re-evaluation of the 1991 Pinatubo eruption using SO 2 derived from HIRS/2 using the Prata fit method and compared it to SO 2 derived from TOMS measurements.They were found to be broadly consistent.The Prata fit method works sufficiently well to suggest that the 7.3 µm SO 2 feature it uses is robust enough to make further exploitation more refined.Use of information arising from other HIRS channels would constitute an improvement to the Prata fit method, as multiple-wavelength information can be used to diagnose attributes of the atmospheric profile such as temperature and the presence of cloud.This problem is well suited to an optimal-estimation (OE) retrieval, which would incorporate a forward model (FM) of sufficient complexity to represent these atmospheric attributes.As with the Prata fit, unavoidably it will require some estimate of the altitude of an SO 2 plume.G. M. Miles et al.: Retrieval of volcanic SO 2 from HIRS/2 using optimal estimation 1.3 Outline of paper In Sect.2, an OE retrieval algorithm methodology to extend the Prata fit method is presented.Section 3 comprises an error study and presents results of retrievals from simulated measurements in order to elucidate the sensitivity of the algorithm and potential sources of error.Section 4 presents a case study of the 1991 Cerro Hudson eruption, where the algorithm is applied to real data and new eruption mass estimates are evaluated and compared to existing mass estimates from other instruments/methods.In Sect. 5 the results are discussed and further work is suggested.

Retrieval algorithm and forward model
The HIRS/2 measurements used here are all-sky brightness temperatures from the instrument aboard NOAA 11.This was selected to demonstrate the capability of this version of the instrument with only one channel that is sensitive to SO 2 and two window channels that have some potential to be used to flag cloud and under some circumstances ash if required (although only one is used here directly).The brightness temperatures are a product derived from the raw voltage measurements via a radiance and brightness temperature conversion and have been subject to calibration factors and some basic quality control.Further information about the instrument is available from NOAA (1981) and elsewhere.The data format contains the time in seconds from midnight of the measurement, the solar zenith angle, 19 IR channel brightness temperatures, one visible channel albedo, latitude, longitude, satellite altitude, line number for each orbit and the scan position (see Table 2).
Retrievals are obtained using the Levenburg-Marquardt minimisation method after Rodgers (2000), and the full optimal-estimation scheme used here is described in detail in Miles et al. (2015).The retrieval uses three HIRS/2 channels to derive three products: the SO 2 column, a scaling factor for a water vapour profile and effective cloud top pressure.The 7.3 µm channel is sensitive to both water vapour and SO 2 .This channel may be said to saturate for SO 2 columns above 600 DU, where significant increases in SO 2 result in small changes in channel BT below the envelope of the channel noise and other error terms.The weighting function for water vapour of the 6.8 µm channel peaks at around 500 hPa (around 5 km) and as such would have some sensitivity to the region where the vast majority of the water vapour in the column resides.To represent both channels accurately, some knowledge of cloud is required, which may be gleaned from the 11.1 µm channel window channel.This channel is highly sensitive to the emitting temperature of the lowest surface it observes (be it cloud or the surface); thus with some knowledge of the surface and atmospheric temperature profile it is possible to obtain an estimate of cloud top height (CTH).Other atmospheric gases not retrieved but that contribute appreciably to channel brightness temperature are represented in the forward model by a climatological value.The potential error that this can introduce is incorporated into the estimate of forward-model error.
Radiative Transfer for TOVS (RTTOV) is a radiative transfer model (developed by the UK Met Office; Saunders et al., 1999;ECMWF, 2001) designed to simulate the instruments of TOVS including HIRS/2, and it is used extensively (particularly for assimilation) because of its speed.It calculates layer transmittances for a variety of trace gas species using look-up tables of parameterised regression coefficients for a range of temperatures and pressures.It has been further developed since the TOVS system was first deployed, and version 10 is used here.RTTOV also has the functionality to compute partial derivatives.
RTTOV estimates channel brightness temperature based on pre-calculated coefficients for layer transmittances that are generated for a range of atmospheric profiles.As such, it is extremely fast, but as it stands it does not incorporate any representation of SO 2 other than at a very low climatological value.To alter the transmittance model to include SO 2 would require substantial re-working of program code.It is possible to calculate a set of predictor coefficients for SO 2 and incorporate them within RTTOV by replacing the properties of another gaseous species that has negligible impact on the total column transmittance within the selected HIRS/2 channels (in this case, carbon monoxide).The coefficients were generated by a "training" methodology using an extensive range of specimen atmospheric profiles, where the SO 2 was represented from very low/background levels to very large perturbations, after Matricardi (2008Matricardi ( , 2010) ) and Siddans (2011).This approach retains the speed and accuracy offered by RTTOV and enables the model to be used to represent atmospheric gases for future instruments not already catered for (ECMWF, 2001).
For this work, the predictors were trained using profiles with up to 300 DU.Some care is required in the generation of these coefficients for SO 2 .They are required to be limited to those that represent a first-order relationship with SO 2 since the more complicated (higher-order) predictors caused erroneous results.This is thought to be a result of both the dynamic range that SO 2 can exhibit in a volcanically perturbed atmosphere and the fact that RTTOV was not explicitly designed to model SO 2 for this instrument.The cost in terms of accuracy over this range of SO 2 is shown to be small, as demonstrated in Fig. 2 and to be discussed in detail later.
The column retrieval developed here uses atmospheric profiles from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim product (Dee et al., 2011) to represent atmospheric properties other than SO 2 , or as a first guess in terms of the water vapour profile.These contain profiles on a pressure grid of 37 levels from 1000 to 1 hPa.RTTOV is capable of generating weighting func- tions, but they refer to the sensitivity of the simulated measurements to perturbations in the atmospheric profile, rather than directly to changes in state vector.As a result, these are evaluated numerically in the forward model by successive FM calls where each element of the state vector is fractionally perturbed in turn.RTTOV has certain physical limits for its input values, and when occasionally the predicted updated state lies outside these they are constrained in the FM by the physical limits that RTTOV will accept, or that are appropriate for the forward model.These are 0.01 to 800 DU for SO 2 , 1 × 10 −6 to 16 times the column water amount predicted by ECMWF, and a maximum cloud top height of 16 km (a conservative upper limit for tropopause height).The weighting functions are allowed to make linear extrapolations beyond these limits, allowing the retrieval more freedom, but unphysical profiles are suppressed with quality control of the derived products (discussed later).

Profile definition in forward model
In the absence of any further information, an effective SO 2 profile must be represented in the forward model.The threeelement state vector comprises a scaling factor for the SO 2 profile, a scaling factor for a water vapour profile and a cloud top pressure.A volcanic SO 2 perturbation is represented by a vertically localised triangular profile.This triangular profile is normalised to have an integrated mass of 1 DU.This was partly done to ease interpretation, since the retrieved scaling factor would be approximately equal to the total amount of SO 2 in the column.The rest of the profile is prescribed by a background SO 2 volume mixing ratio climatology, the total column mass of which is less than 1 DU.In the forward model, a scaling factor applies to a specified height region of the SO 2 profile, scaling all elements within and none out-side this.The expected region of the volcanic plume is estimated using ancillary information, such as lidar or results from modelling of the eruption available in the literature.The retrieval sensitivity to how well thickness and altitude of the plume are modelled as compared to the true state is evaluated using retrievals from simulated measurements.These are detailed in Sect.3.
In an analogous way to SO 2 , H 2 O is represented in the state vector by a profile scaling factor, but it applies to the entire profile rather than a localised height region.The profiles used for retrieval are those collocated from the ECMWF ERA-Interim product for a given HIRS/2 pixel (which represents the best guess for the state), but in principle any climatological profile can be used.In the case where a scaling factor is close to 1, it would indicate that the H 2 O profile is similar to that which produced the measurement.
The third element of the state vector is CTH, or specifically the geopotential height at an equivalent pressure level.It was found that the speed of convergence was significantly reduced if the initial guess of cloud top pressure was reasonably accurate.As such, this is derived before the retrieval using interpolation between calls to a radiative transfer model that simulates the 11.1 µm channel BT for 0-10 km (using associated ECMWF ERA-Interim temperature profile), and it includes a test for temperature inversions.

Error
An estimate of forward-model error was calculated using the Reference Forward Model (RFM) -a line-by-line radiative transfer model (Dudhia, 2002), discussed further in Sect.3. The estimate accounts for inaccuracies that arise due to modelling the atmosphere at reduced spectral resolution, limited vertical resolution (100 m versus 1 km as used in the forward model outside the region of the SO 2 perturbation), inclusion of non-retrieved trace gases at a climatological level or their preclusion entirely, relative to a reference case.This yields a channel quantity (in brightness temperature) that is combined in quadrature with the noise-equivalent differential radiance for each channel and is thus incorporated into measurement noise for the purposes of the retrieval.The a priori error associated with cloud height is 10 km.The a priori error for water vapour is based on the variance of water vapour in the ECMWF atmospheric training profiles discussed above relative to the mean.

Estimation of SO 2 and H 2 O covariance for HIRS/2
Establishing an appropriate SO 2 a priori error is potentially a non-trivial issue with regard to a retrieval problem where the measurements have relatively little sensitivity.A volcanically perturbed SO 2 profile can contain 2 or 3 orders in magnitude more than a background profile, and at the centre of a large plume this can be even more.A good a priori error gives the retrieval the freedom to find a correct minimum in cost www.atmos-meas-tech.net/10/2687/2017/space and can restrict it from converging on a solution that is unphysical.The variance for a background profile would be very small, as opposed to a profile where SO 2 is expected, which would be very large.If there is sufficient information contained within the measurements, one would conventionally use a variance that spans both scenarios.This results in a poor constraint for an ill-posed problem but is necessarily used here, where a first guess/a priori error of 100 DU is used and a prior variance is the first guess squared.A column amount of 100 DU may represent that from a large, explosive volcanic eruption.Pinatubo, for example, yielded column amounts of 350-500 DU (depending upon instrument) after 24 h, which reduced to 100 DU after 7 days (Carn et al., 2005).The OMI instrument (see Table 1) captured column amounts of around 200 DU after the 2008 eruption of Kasatochi (Prata et al., 2010).
Early results of the retrieval scheme run with real measurements revealed that there were many "false positives" of SO 2 retrieved.Their structure indicated that they were related to the presence of water vapour, or errors in the fit for water vapour.This indicated the degree of covariance between SO 2 and water vapour which had to be incorporated into the retrieval since the 7.3 µm channel is sensitive to both water vapour and SO 2 .
The retrieval was applied to one day of "clean" measurements in the Southern Hemisphere, where no volcanically perturbed profiles were expected.The retrieval was forced not to retrieve SO 2 by artificially constraining the a priori variance, but nonetheless small amounts of SO 2 are retrieved from that channel because of inadequacies in characterising the water vapour.The brightness temperature fit residuals in the SO 2 channel were very small, but it is expected that nearly all of the SO 2 being retrieved on this day is being falsely attributed.The standard deviation of the 7.3 µm channel brightness temperatures fit residual in the retrieval of 0.92 K constitutes an estimate of the "real-world" error covariance of water vapour with SO 2 for this instrument.This is incorporated by adding it in quadrature to the forward-model error for this channel and resulted in a significant reduction in the occurrence of false positives.

Error study: retrievals from simulated measurements
There are some sources of error that can be incorporated into and dealt with by the retrieval.These include measurement noise, the presence of cloud or ash, SO 2 / H 2 O covariance and an estimate of forward-model error discussed above.
The main sources of error that cannot be adequately represented in the forward model are errors that impact ill-posed nadir SO 2 column retrievals in general.These are incorrect height assignment of the SO 2 plume; incorrect thickness in the plume represented in the forward model; and, particularly in the case of infrared measurements, sensitivity to the presence of cloud and/or water vapour.Their relative impacts vary, and the sensitivity of the solution to them can be quantified using simulations.It should be noted that some of these errors (plume height and profile shape) cannot often be known at the time of retrieval, and as such the actual impact on the retrieval result also cannot be known.They are investigated here in order to give a general indication as to the potential error that can be associated with the results, to give a window of confidence.Others, such as the impact of cloud or ash on the retrieved SO 2 error, can be investigated for use in quality control.

Spectral precision of forward model
In order to assess the accuracy of the RTTOV-based fast column retrieval forward model, it is compared to simulations from a model with a higher accuracy.The RFM is a lineby-line radiative transfer model (Dudhia, 2002) capable of modelling the atmosphere at a spectral resolution of up to 0.0001 cm −1 .The RFM is not suitable for the forward model because it is computationally expensive and does not inherently represent any effects of cloud or ash. Figure 2 shows the results of column retrievals from HIRS/2 channel BTs simulated by the RFM, using a sample ERA-Interim cloud-free meteorology (temperature and water vapour profiles) at 0 and 60 • S latitude and 0 • W longitude, where only the column amount of SO 2 is changed in the simulation.It also shows the SO 2 fit by the Prata fit method.The Prata fit method does not fit SO 2 below 5 DU, which depending upon the atmospheric state can be equivalent to an observed brightness temperature difference of up to 4 K.The bias of the Prata fit has a dependence upon latitude, primarily because of the different amount of water vapour in the profile at the two latitudes shown here.The column retrieval has a very small bias that only becomes perceptible at SO 2 loadings approaching 200 DU, at which point it is of the order of < 5 DU.

Sensitivity to forward-model representation of SO 2 plume
Both the altitude and amount of SO 2 affect the 7.3 µm channel brightness temperature, but as there is only one channel sensitive to SO 2 on NOAA 11 considered here, there is at most one piece of information that can be retrieved for SO 2 .Therefore, for an accurate retrieval of SO 2 column, it is important to have some knowledge of the plume altitude or its vertical profile.The column retrieval developed here requires some information of the height of the SO 2 , but this can be subject to uncertainty and may change with time.As such, the sensitivity of the retrieval to errors associated with plume height and specification must be examined.

Altitude
Measurements were simulated for a plume at a range of altitudes from 8 to 18 km.Figure 3 shows the impact on Figure 3.A measurement was simulated for a volcanic plume of triangular profile centred at a range of altitudes, for a range of total column amounts.A retrieval is then performed where the plume is assumed to be at 12 km.The fractional difference, or error, is plotted.
the retrieved SO 2 column at a specified, fixed altitude of 12 km as a fraction of the true column at these altitudes.
Errors range from typically ±0 to 30 % for most column amounts up to 100 DU and increase for larger amounts and for particular altitudes.Whilst the specific error may be statedependent (upon meteorological conditions, specifically the water vapour profile), these simulations do give a general indication as to the magnitude of error that can result from incorrect height assignment of the volcanic plume in the forward model.This is the largest source of error in the OE column retrieval (and the Prata fit method) and is made more challenging because there is a dependency of the error on column amount.Since height assignment errors cannot be known, such simulations can at least give a general indication of potential uncertainty of retrieved amounts, depending on the quality of information available regarding altitude of volcanic SO 2 .It is clear therefore that good prior knowledge of the SO 2 plume altitude is necessary for accurate retrieval or fit of SO 2 column amounts from HIRS/2.The performance of the column fit was also directly assessed against a line-by-line model (RFM) for plume altitudes from 8 to 18 km (where the plume height assignment used in the retrieval was the same as that used in the measurement simulated by the RFM), and it was found that for altitudes of over 17 km the column fit was unable to retrieve SO 2 columns less than 30 DU, but in all other cases true clear-sky column amounts were retrieved accurately from simulated measurements.

Profile shape and plume thickness
Figure 4 shows the consequences that can result from retrieving the volcanic plume with a fixed profile shape that repre- sents the thickness of the plume incorrectly.Measurements were simulated using a triangular profile centred at 12 km but with baselines of 1 and 4 km.They were then used in the retrieval with a fixed profile shape with a triangular perturbation also centred at 12 km but with a baseline of 2 km (thought to be the best representation of the plume used in the case study in Sect.4).The retrieval simulations suggest that errors are larger when the plume thickness is overestimated (typically 13 %), with only small inaccuracies introduced when the plume thickness is underestimated (less than 2 %).The modelled cloud top height was 3 km in all cases.It is therefore possible that an underestimate of plume thickness would result in smaller errors.

Sensitivity of retrieval scheme to cloud and ash
Some understanding must be obtained of how the column retrieval forward model behaves in the presence of ash and cloud of different type.The forward model fits a cloud top pressure using the 11.1 µm channel, which is expected to work well for most scenes with cloud in the troposphere.The effect of cloud on the other channels is examined here using a cloud model, the Oxford-RAL Retrieval of Aerosol and Cloud (ORAC) model.The model is described in detail by Poulsen et al. (2012), where it was used as part of an optimalestimation retrieval of cloud properties for the Along-Track Scanning Radiometer (ATSR) by simulating radiances in a combination of visible, near infrared (NIR) and IR channels.The model parameterises a cloudy scene by ascribing cloud phase, effective radius of a size distribution, the 0.55 µm opwww.atmos-meas-tech.net/10/2687/2017/tical depth and a cloud top pressure.It uses the plane-parallel approximation and models cloud as a single layer.The model represents trace gases at a background climatological level.The system can also be used to retrieve ash plume properties: plume height, optical thickness and ash particle effective radius (McGarragh et al., 2017).
HIRS/2 measurements were simulated for a range of liquid and ice cloud and ash optical depths, for a range of effective radii and at a range of altitudes when no volcanic SO 2 was present.These channel brightness temperatures were then used to retrieve SO 2 to identify where this resulted in an erroneous fit.
An example is shown in Fig. 5, which shows that for liquid water clouds above 5 km the column retrieval erroneously retrieves some SO 2 when there is none, the water vapour and cloud top height become inaccurate and the fit cost begins to increase.The results indicated that low optical depth or effective radii for cloud or aerosol can result in poor fitting of the measurements, resulting in both an underestimate of cloud top pressure with false positives of SO 2 and an overestimation of water vapour.This yields a crucial quality control threshold where retrieved cloud top altitudes of greater than 5-6 km should not be trusted, as they are likely to result in spurious detection of SO 2 and a high retrieval cost.This may imply that very thin cloud beneath 5 km (or incorrectly retrieved to be) could still contribute to poor fitting of the measurements.

Quality control
The results of the column retrieval must be subject to some quality control.In addition to the disregard of non-converged and converged pixels with cloud retrieved at an altitude greater than 5 km, a retrieved column is only considered useful if the error is less than the retrieved amount.Quality control becomes very important when erupted plumes are used to calculate total erupted mass, where even a small amount of noise can yield a biased mass total.For the purposes of gridding or summing pixels for deriving a global/plume mass estimate, a minimum retrieved SO 2 threshold may be applied in deference to the lower detection limit of the retrieval, in order to avoid spurious low values that the retrieval should not be sensitive to, such as those relating to water vapour or cloud that are not accounted for in either the error covariance or the forward model.An effective way of obtaining this quantitatively is to apply a 2 or 3σ test, where sigma is the standard deviation of the retrieved SO 2 on a day when no volcanic SO 2 is expected to be present.This threshold gives statistical confidence that a value above it is significantly distinct from the noise above the 95th or 97th percentile.The sigma threshold for 6 August 1991 (a day when there was no SO 2 present in the region relating to the case study in Sect.4) was 2.7 DU and is probably a lower estimate of the detection limit of the HIRS/2 SO 2 column retrieval at the mid-latitudes.Multiples The eruption was estimated to be 10-20 times smaller than Pinatubo in terms of SO 2 that was expected to be emitted.In this sense, as well as being a non-equatorial eruption, it has similarities to the 2008 Kasatochi eruption in the Northern Hemisphere.It is selected here as a case study because it was a relatively large eruption that has not been studied exhaustively, and it is a very good example of an eruption in recent satellite history which only TOMS observed with any significance that can benefit from application of this technique.
At the time of the 1991 eruption, the only satellite available that could detect SO 2 with any demonstrated accuracy was TOMS.The Microwave Limb Sounder, a contemporaneous instrument that observed SO 2 from Pinatubo at a higher altitude, produced noisy results in the lower stratosphere at this latitude (Read et al., 1993).In addition, contemporary lidar measurements of the Hudson plume were made at the CSIRO (Commonwealth Scientific and Industrial Research Organisation) Division of Atmospheric Research, at Melbourne, Australia (38 • S, 145 • E) (Young et al., 1992;Barton et al., 1992).These measurements are sensitive to ash, sulfate aerosol and meteorological (water) cloud.The backscatter profiles tend to indicate peaks at around and above 20 km, and frequently at 10-13 km.The higher peak is attributed to aerosol from the Pinatubo eruption.Young et al. (1992) interpret the majority of observations that are thought to include Hudson material as the feature at 12 km in October, with variable cirrus at 10 km.It is reported by the authors that the plume was observed consistently from 28 August until December 1991 between 10 and 13 km, with a decreasing scattering ratio.The relative proportions that contribute to the backscatter measured are expected to be dominated by ash in the first few weeks after the eruption.Little ash is expected to be present after a month beyond the eruption, but by this time the vast majority of the SO 2 will have oxidised into aerosol.Whilst lidar is not sensitive to the presence of gaseous SO 2 , inferences can be drawn from the height of the aerosol it eventually becomes.In this case the lidar information is considered to be a valuable starting point as a guide for estimating the cloud height of the SO 2 , in the context of other information.As well as some ground observations, the Hudson eruption was sensed remotely by AVHRR (ash), lidar (sulfate aerosol) and incidentally by an aircraft (Barton et al., 1992).Hofmann et al. (1992) reported possible exacerbation of Antarctic ozone depletion of 10-20 % of total column due to the presence of Hudson aerosol in the lower stratosphere for September 1991.The anomalous depletion occurred within the polar vortex predominantly at 11-13 and 25-30 km, the respective altitudes of the Hudson and Pinatubo aerosols.
The transport of the Hudson volcanic plume was first numerically modelled by Barton et al. (1992), to reasonably good agreement with satellite and lidar observations.The plume was also modelled using an isentropic trajectory model, initiated by TOMS observations of SO 2 (Schoeberl et al., 1993).These models showed good spatial agreement with observations for the first 8 days after the eruption, which is an indication that the height assignment of the erupted plume was accurate within the models.The most explosive eruption began and ended on 15 August.It was at this stage of its eruptive phase that the majority of the material was injected into the stratosphere (Constantine et al., 2000).

Results
Using all of this information, the Hudson plume is modelled as a triangular peaked profile with a baseline of 2 km between 11 and 13 km, peaking at 12 km.Figure 6 shows an example of the SO 2 retrieval applied to a day of data on 15 August 1991 and its associated retrieval error.Figure 7 shows results for the same day as Fig. 6, but for the other elements of the state vector: the retrieved water vapour scaling factor and cloud top height (with their associated retrieved errors).Only high-cost and convergence criteria have been applied.In general, the retrieved values of cloud top height have very small errors.For the water vapour scaling factor, the largest errors occur in the presence of high or thick cloud, which is expected.As shown in Sect.3, the cloud model simulations suggested that the retrieval struggles in the presence of high cloud and can on occasion fit spuriously enhanced SO 2 , potentially because it results in a poor estimate of water vapour in the correspondingly colder scene.Regions of very high water vapour scaling factor result in very high errors in retrieved SO 2 , and data with cloud top height greater than 5 km are not considered reliable for SO 2 .
Figure 8 shows 9 days of retrieved SO 2 from the 1991 Cerro Hudson eruption following the largest eruption phase on 15 August.The eruption began on 8 August, emitting smaller amounts of SO 2 into the upper troposphere-lower stratosphere, which can be seen as already present in the path of the main plume on subsequent days.The multiple sampling of the plume by successive orbits (day and night) is quite apparent, particularly as the plume becomes more distorted after 20 August.

Plume mass estimate
The simplest method to estimate the total erupted mass or mass present in a volcanic plume is to take the sum of the representative footprint areas of the satellite that measured SO 2 .This method presents several problems relating to sampling of a volcanic plume; particularly with an infrared instrument that measures both night and day and that could sample the plume more than once, orbits may partially sample the plume in any one swath, and the plume will move constantly between sampling.Alternatively, gridding averages the data into grid boxes on a latitude-longitude grid.Some care must be taken to account for whether or not the gridded data are representative of the data resolution, and keeping track of bins with no data can be a way to estimate undersampling.Guo et al. (2004) used two methods of gridding data, that of kriging for TOMS data and nearest-neighbour interpolation for HIRS/2 (Prata fit method) to account for larger spatial gaps between points.These methods either impose statistical methods or manually introduce information based on assumptions.Whilst both can be utilised in such a way as to indicate an estimate of the error or uncertainty that this introduces, mass estimates presented here are only based on the sum of equivalent contiguous footprints represented by each HIRS ellipse.
Furthermore, if gridding is used, in order to ensure that the data are sampled fairly, the orbits should first be split into ascending and descending nodes, with care taken regarding where a plume is in relation to the date line.This is in an effort to minimise recording the same data point twice when the plume has moved by the time the region is sampled again.Other methods are available but often require a model or further ancillary information.

Comparative measurements of SO 2
The plume mass estimate for the HIRS/2 SO 2 retrievals for the Cerro Hudson eruption may be qualitatively compared to the figures for TOMS within Constantine et al. (2000).Total erupted mass estimates given can be directly compared, as shown in Table 3, although the methodology by which the estimates were derived differs.Spatially, HIRS/2 has the advantage of a smaller footprint than that of TOMS (instrument field of view (IFOV): 1.25 • × 1.25 • / 17.4 km × 17.4 km versus 3 • × 3 • / 50 km × 50 km), but the TOMS swath is 50 % wider (3000 km).For a case such as the Hudson plume, TOMS is more likely to capture the entire plume in one orbit swath and sample it only once, which on the one hand  greatly reduces ambiguity in deriving total plume mass, but on the other hand the frequency of observation is reduced, and sometimes only part of the plume is captured.As reported by Constantine et al. (2000), this was sometimes the case, and a "best" estimate of the TOMS data was used to contribute to the values in Table 3.The erupted mass estimates given in Table 3 that relate to HIRS/2 are the sum of equivalent footprint areas, from nodes that capture the most of the SO 2 plume present each day.Figures are rounded to reflect probable accuracy.For the total eruptive period, this method has yielded a total erupted SO 2 mass estimate of 2300 kT with an averaged retrieved error of 27 %.This error does not incorporate error that arises from uncertainty in the height of the SO 2 in the forwardmodelled plume (as demonstrated in Sect.3), or error that might arise from discounting pixels where SO 2 was retrieved below the 3σ threshold.It does not account for absent scan lines due to instrument calibration, so it should be considered a lower limit.As previously discussed, a good estimate of plume height is an unavoidable requirement in SO 2 detection with an instrument with only one channel sensitive to atmospheric SO 2 .In the case of this work, height assignment error of ±1 km introduces a mass-dependent bias of between 5 and 20 % for a given pixel depending upon where in the atmosphere the plume is located.For TOMS, the approximate error suggested for the total erupted mass estimate is 30 % (Krueger et al., 1995;Constantine et al., 2000).
The TOMS algorithms used in Constantine et al. ( 2000) have been recently updated, and a brief comparison is presented here to some initial data from an updated TOMS algorithm.This algorithm exploits the way ozone and sulfur dioxide both strongly absorb UV radiation.The new TOMS algorithm builds on the early heritage of BUV algorithms (Krueger et al., 1995).These algorithms retrieve both O 3 and SO 2 by taking advantage of the large SO 2 / O 3 cross section ratio (CRS) differences in the gas-absorbing bands.This approach constructs radiance tables using a forward model that accounts for both the O 3 and SO 2 cross sections.The new algorithm uses the 317 nm channel to retrieve SO 2 (CRS ∼ 2.5); the 331 nm channel to retrieve O 3 (CRS ∼ 0.15); and the channel at 340 nm to retrieve the spectral dependence, dR / dλ.This methodology further applies a small second-order correction that accounts for nonorthogonality between the SO 2 and O 3 channels.
A 1-week composite of retrieved SO 2 for both instruments is shown in Fig. 9, where SO 2 from the main eruptive phase can be seen circumnavigating the hemisphere.There is clear complementarity between the instruments in terms of absolute amount retrieved and characterisation of the plume.The smaller pixel size of HIRS and more frequent sampling enable the plume to be observed in finer detail; however the wider swath of TOMS frequently captures more of the plume in one swath.
For a more detailed comparison, two orbits during the 1991 Hudson eruption are considered where the plume is al-  most fully sampled by both instruments, as shown in Fig. 10.The pixels in the region of the plume were also relatively cloud-free or had low cloud during the observation.The geographical bounds considered for the mass estimate are between −53 and −45 • in latitude and 10 to 60 • in longitude.Using the method of summing over mass and area discussed previously, the mass of the plume represented here by HIRS/2 and TOMS is calculated to be 1398 and 1540 kT respectively, after quality control has been applied.The missing four scan lines due to a HIRS calibration phase that coincide with the plume in the region of high concentration sug-gest the HIRS estimate is an underestimate.It is apparent that HIRS/2 is potentially more sensitive to lower amounts of SO 2 .It is challenging to directly compare the SO 2 retrieved by two instruments with differing footprint sizes.Gridding might offer an alternative method of plume mass estimate, but selection of the most appropriate grid box size relative to the pixels of each instrument coupled with the small size of the plume with a strong SO 2 concentration gradient make it a challenge for such a comparison to be equitable and account for instrument attributes.A comparison involving gridding for a larger eruption (c.f.Pinatubo) would be less problematic.

E-folding time
The e-folding time for erupted SO 2 is a measure of the residency of the material in the atmosphere and is affected by the height the material reaches and, in the case of very large eruptions, the amount itself.It is also affected by wind shear (horizontal and vertical) and humidity, which affects the rate at which the SO 2 is oxidised and sulfate aerosols grow.The measure is more suited to large eruptions (e.g.El Chichn in 1982 or Pinatubo in 1991), in terms of inferring effects upon radiative forcing, about which Miles et al. (2004) and other works are concerned.This is because the amount and height that such eruptions reach in the stratosphere give the SO 2 sufficient time to become globally mixed and as such affect the radiative forcing globally.Equation (1) describes the process of exponential decay, where N (t) is a quantity at time t, N 0 is the initial quantity at time t = 0 and λ is the decay constant.
The e-folding time, the time in which the initial quantity is reduced to 1/e of its initial value, is given by the reciprocal of the decay constant.Using approximate values from the mass estimates derived from Fig. 9 where the total SO 2 can be said to drop from around 1500 kT (the total mass present on 17 August 1991 associated with main plume) to 500 kT 18 days later, this yields an e-folding time of around 16 days.Two days after the largest plume was erupted is used here to minimise potential obscuration of the plume by the coincident presence of thick ash.In reality the total mass observed does not decay smoothly but has noise due to the fact that the plume is not always perfectly sampled, and the number of retrieved pixels is excluded due to the presence of high or thick cloud or ash varies.The variability of the mass estimates and the associated retrieval error make only an estimate appropriate for this approach, but it is not considered to be an unreasonable one.If the e-folding time is calculated for the extremes of the retrieved error bounds of the mass estimates, the e-folding time is 10 days at a minimum and 35 days at its shallowest descent, but these are considered to be overly generous bounds by this method.This case is complicated by the fact that about 30 % of the SO 2 released by Hudson was erupted over the 7 days before the main eruption on 15 August, making the calculation of the decay subject to further uncertainty.The e-folding time for this SO 2 plume as estimated by Constantine et al. ( 2000) is around 15 days, but they state that this is algorithm-dependent.These estimates are somewhat smaller than the e-folding times for the larger eruptions (e.g.Pinatubo), which is to be expected due to the considerably lower altitude of the Hudson plume.More recently, Carn et al. (2016) estimated the e-folding time of Cerro Hudson to be ∼ 7 days, based on mass estimates from TOMS (Constantine et al., 2000).They attribute this anomalously short e-folding time to the late Southern Hemisphere winter timing of the eruption.However, since Constantine et al. (2000) estimate nearly twice the initial total mass (4000 kT) than that observed by HIRS/2 in this work (and the subsequent TOMS algorithm discussed here), it is possible that the inconsistency in e-folding times could be due to an overestimate of initial erupted mass from the original TOMS algorithms.Total mass estimates (and therefore efolding time estimate) would be improved greatly in accuracy if the HIRS/2 instruments aboard NOAA 10 and NOAA 12 that were also present were used to result in very comprehensive sampling of this eruption.

Discussion
This OE column retrieval finds a new total erupted mass estimate for the 1991 eruption of Cerro Hudson of 2300 ± 600 kT from the HIRS/2 instrument aboard NOAA 11, where the error is the retrieval noise.This does not incorporate any error from plume altitude estimation, but the potential impact has been quantified by forwardmodel simulations.This total mass estimate is lower than that of TOMS (Constantine et al., 2000) and that of Carn et al. (2016) but higher than that derived in a similar way using the methodology of Prata et al. (2003) for HIRS/2.Reasons for this include (but are not limited to) differences in sampling, height sensitivity, instrument differences and attributes or accuracies of the forward model or fit employed in SO 2 detection.From the comparison with the new TOMS algorithm, the HIRS/2 results presented here are highly consistent, and further quantitative comparison, for this eruption in particular, is desirable.The retrieval precision demonstrated in this case study is slightly smaller (∼ 3 DU) than that proposed for the TOMS instrument (6-7 DU).As such, with the increased sampling of the IR instrument it is apparent that HIRS/2 can offer a positive contribution to the atmospheric SO 2 emission record from explosive volcanic eruptions up to and beyond the launch of GOME and other satellites that followed.Moreover, benefits of the optimal-estimation approach over and above the more rapid but limited brightness temperature difference method are significant.They include a quantified error on individual pixel retrieved values, latitudinal varia-tion in accuracy, diagnostic indicators of the retrieval performance and goodness of fit and treatment of cloud and water vapour consistent to the retrieval of SO 2 .When summing mass over a large number of pixels, the precision that these afford becomes increasingly important.Issues that remain are those endemic to ill-posed problems where there is only one piece of information on SO 2 available and only limited information about the height or shape of the profile of a volcanic plume.It is conceivable that further progress might be made by using HIRS/2 aboard NOAA 10 and 12 with the addition of the 8.6 µm channel in ash-free pixels.
There are clear opportunities for extending this work.In particular, as the HIRS/2 instrument was present aboard a number of the NOAA platform series and often simultaneously flown (NOAA 10, 11 and 12 were all in orbit at the time of the Cerro Hudson eruption), there is the possibility to fully characterise eruptions with very high temporal sampling.More rigorous methods for interpolation, sampling and gridding the data can also be used to reduce errors in the total mass estimates.The application of further tools such as chemistry transport or trajectory models for understanding plume evolution would be better constrained by the availability of more measurements.
The first HIRS instrument was flown aboard TIrOS-N in 1978, and there are almost continual data available to the present, and for the foreseeable future of the Met-Op series of satellites, enabling a potential data set spanning 40+ years.Generating an SO 2 data set for the duration would be an opportunity to maximise the value and legacy of the satellite data.Such a data set, with an accompanying error covariance estimate, could be used as input to a climate model to better assess the effects of large volcanic eruptions on the radiative balance of the atmosphere.For much of the latter half of that period, there are (and will be) other satellite instruments capable of measuring SO 2 in the limb and the nadir, in particular high-resolution spectrometers with very much enhanced accuracy and precision that will provide correlative information about the quality of the HIRS/2 SO 2 column retrievals that may be considered in retrospective terms.There is also a break in the TOMS record during 1995-1996 that can be filled by HIRS/2 estimates.
It would be highly desirable to extend comparisons from this eruption with TOMS SO 2 in general, comparing a longer record by both instruments for other eruptions, since both provide a unique record of SO 2 potentially spanning many decades.Satellite records of this length for climatologically important trace gases are rare and would also provide further constraint to volcanic SO 2 emissions in coupled chemistryclimate models.

Figure 1 .
Figure1.Transmission spectra of H 2 O and SO 2 simulated from Southern Hemisphere mid-latitude water ECMWF ERA-Interim background vapour profile using the RFM (see text).The SO 2 spectra were simulated using triangular profiles to represent column amounts of 1 and 300 DU, as used in the forward model.

Figure 2 .
Figure 2. Retrievals based on simulations by a line-by-line model (RFM), with synthetic measurement noise.The error bars for the column retrieval are the retrieved errors.These simulations use temperature and water vapour from a cloud-free ECMWF ERA-Interim atmosphere on 15 August 1991, for a grid box centred at 0 • N and −60 • N latitude and 0 • E longitude.The vertical bars show the retrieved error for the column retrieval.No error estimates are possible for the Prata fit method.

Figure 4 .
Figure 4.The black line indicates how columns from 0.1 to 200 DUare retrieved on a fixed grid with a scalable triangular profile with base, mid-point and top at 11, 12 and 13 km respectively, when the true profile shape is given by a triangular profile at 11.5, 12 and 12.5 km, effectively overestimating the thickness of the plume.The red line shows the equivalent result for an underestimate of the plume thickness, the real profile given by 10, 12 and 14 km.The dotted lines show the bounds of retrieved error in each case.The dashed line is x = y, shown for clarity.

Figure 5 .
Figure 5.The top left plot shows retrieved cloud top height as a function of "true" cloud top height as simulated by the cloud model.Black symbols indicate that the retrieval converged, and purple indicates that it did not.The top right plot is of the fit residual (measurement minus fit) in the 11.1 µm channel.The bottom left plot shows the retrieved SO 2 as a function of the cloud top height in the cloud model, and the bottom right the equivalent for the water vapour scaling factor.

Figure 6 .Figure 7 .
Figure 6.Retrieved SO 2 columns for 15 August 1991, and retrieved error for orbits that day.Erupted SO 2 from the start of the eruptive phase (from 8 August 1991) is evident ahead of the larger plume emitted on 15 August.Data are screened at the 2σ level (5.4 DU).

Figure 8 .
Figure8.Progression of main erupted plume from 15 August 1991, using all orbits (day and night) from HIRS/2 NOAA 11.The eruption began with smaller amounts emitted from 8 August, which are apparent on the 15th and disassociated from the main plume.The plume's transport between observations is evident, particularly from 21 August, where it is captured multiple times by multiple swaths.Data have been screened at the 3σ level (8.1 DU) for clarity of the main plume.

Figure 9 .
Figure 9. Seven-day composite of retrieved SO 2 from 15 to 21 August 1991.For clarity in comparison, TOMS data are screened to have a minimum value of 15 DU, and HIRS/2 data use 3σ (7.1 DU).

Figure 10 .
Figure10.The main Hudson plume on 17 August 1991 as observed in orbits 5 and 6 by HIRS/2 and 64695 and 64696 by TOMS, 2 days after the main paroxysmal eruption that occurred on 15 August.Four scan lines in the HIRS/2 panels are missing due to a routine calibration phase in which no data are provided.HIRS and TOMS data are both screened at the quality level of 2σ level (5.4.and 15 DU respectively).

Table 3 .
Total erupted SO 2 rounded estimates for Cerro Hudson.