SO2 Layer Height retrieval from Sentinel-5 Precursor/TROPOMI using FP_ILM

Accurate determination of the location, height and loading of SO2 plumes emitted by volcanic eruptions is essential for aviation safety. The SO2 layer height is furthermore one of the most critical parameters that determine the impact on the climate. Retrievals of SO2 plume height have been carried out using satellite UV backscatter measurements, but until now, such algorithms are very time-consuming. We have developed an extremely fast yet accurate SO2 layer height retrieval using the Full-Physics Inverse Learning Machine (FP_ILM) algorithm. This is the first time the algorithm has been applied to 5 measurements from the TROPOMI instrument on board the Sentinel-5 Precursor platform. In this paper, we demonstrate the ability of the FP_ILM algorithm to retrieve SO2 plume layer heights in near-real-time applications with an accuracy of better than 2 km for SO2 total columns larger than 20 DU. We present SO2 layer height results for a selection of recent volcanic eruptions observed by TROPOMI. 10 Copyright statement.


Introduction
Global satellite observations allow for the timely detection and monitoring of SO 2 emitted from volcanic eruptions, even in remote regions, where no ground-based instruments are installed (see e.g. Fioletov et al. 2013). Satellite measurements of UV earthshine spectra in the wavelength range between 305 and 335 nm provide the highest sensitivity to SO 2 in the Earth's atmosphere. Volcanic eruptions can inject large amounts of sulphur dioxide (SO 2 ) into the atmosphere, where it is either subject 15 to dry and wet deposition within a few days in the troposphere (see e.g. Lee et al. 2011or Myles et al. 2011, or oxidization over a period of several weeks to sulfate aerosols in the stratosphere (see e.g. Robock, 2000;Forster et al., 2007 andvon Glasow et al. 2009). Sulfate aerosols can affect the Earth's radiative forcing and have an impact on clouds (see e.g. McCormick et al., 1995;Robock, 2000 andMalavelle et al. 2017).
Based on UV earthshine measurements, SO 2 vertical column densities (VCDs) can be retrieved easily using for example 20 the Differential Optical Absorption Spectroscopy (DOAS) algorithm, see e.g. Rix et al. (2012), or the Principal Component The determination of the vertical distribution of SO 2 based on satellite data is however challenging, since the information about the vertical distribution is not easy to extract from the spectral signature (see Yang et al. 2009 andNowlan et al. 2011).
The SO 2 loading (VCD) has a direct effect on the optical depth, whereas the layer height (LH) has an indirect effect on the optical depth since it influences both the number of photons passing through the SO 2 layer as well as the layer optical depth due to the temperature dependency of the SO 2 absorption crosssections (Yang et al., 2009). 15 In the operational Sentinel-5 Precursor/TROPOMI product (see Theys et al. 2017), the SO 2 VCDs are hence provided for a set of apriori SO 2 distributions, along with an averaging kernel (AK), in order that the user can calculate the VCD for an arbitrary SO 2 vertical distribution.
To date, SO 2 layer-height retrievals have used computationally demanding direct fitting inversion methods, which are not suitable for NRT applications. For the retrievals based on satellite UV measurements, Yang et al. (2009) and Yang et al. (2010) 20 developed an Extensive Iterative Spectral Fitting (EISF) algorithm for OMI, and Nowlan et al. (2011) introduced an optimal estimation (OE) scheme for GOME-2. For strong volcanic eruptions, the accuracy of the retrieved SO 2 LH using this approach is in the range 0.5-1 km, whereas for small SO 2 absorption it is around 2 km (see Nowlan et al. 2011).
Satellite infrared sounders also offer the opportunity to measure both SO 2 VCDs and LH (see Clarisse et al., 2008;Carboni et al., 2012;Clarisse et al., 2014 andCarboni et al. 2016) using optimal estimation algorithms for IASI. While infrared sounders 25 have a weaker sensitivity to the lower tropospheric SO 2 than UV instruments, the layer height retrievals tend to have a better accuracy, and perform well even for low column amounts (up to the DU-level), however at a low horizontal resolution of about 12 km.
In contrast, the TROPOMI instrument, on board the Sentinel-5 Precursor satellite (S5P) launched on 13 October 2017, has a much higher spatial resolution of 7 × 3.5 km 2 , with the possibility of operation at an even smaller resolution of 5.5 x 3.5 km 2 . 30 This allows us to observe and study SO 2 plumes at an unprecendented level of detail. Data turnover from TROPOMI is very large, and this consideration will require the development of new retrieval schemes for the fast and accurate retrieval of SO 2 layer heights in an operational environment.
To this end, we have developed an algorithm called 'Full-Physics Inverse Learning Machine' (hereafter referred to as FP_ILM) for the retrieval of the SO 2 LH based on satellite UV earthshine spectra. The FP_ILM algorithm has been used 35 2 Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2019-13 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 25 February 2019 c Author(s) 2019. CC BY 4.0 License.
for the retrieval of ozone profile shapes (Xu et al., 2017), the retrieval of surface properties accounting for the BRDF effects (Loyola et al., 2019), and the retrieval of SO 2 LH from GOME-2 (Efremenko et al., 2017). The algorithm creates a mapping between the spectral radiance and SO 2 LH using machine learning methods. The time-consuming training phase of the algorithm using radiative transfer model calculations is performed off-line, and only the inversion operator has to be applied to satellite measurements -this makes the algorithm extremely fast and it can thus be used in near-real time processing environments. In 5 this second paper on the FP_ILM SO 2 LH, we describe some improvements to the algorithm, and we apply it to a number of volcanic eruptions observed by TROPOMI during the commisioning phase (started 7 November 2017) as well as the first operational phase (started April 2018).
The paper is organized as follows: In Section 2 we describe the improved FP_ILM SO 2 LH algorithm. The sensitivity of retrieved SO 2 LHs to a number of different parameters is discussed in Section 3. In Section 4, the FP_ILM is applied to S5P 10 data to retrieve SO 2 LHs for selected volcanic eruptions. We summarize the paper in Section 5.

FP_ILM algorithm
Conceptually, the FP_ILM consists of a training phase, in which the inversion operator is obtained using synthetic data generated with an appropriate radiative transfer (RT) model, and an operational phase, in which the inversion operator is applied to real satellite measurements. The main advantage of the FP_ILM over classical direct fitting approaches is that the 15 time-consuming training phase involving complex RT modeling is performed offline; the inverse operator itself is robust and computationally simple and therefore extremely fast.
During the training phase, the LInearized Discrete Ordinate Radiative Transfer model (LIDORT) with inelastic rotational Raman scattering (RRS) implementation (Spurr et al., 2008) is deployed to compute simulated reflectance spectra in the wavelength range 310 − 335 nm. These spectra depend upon the following n = 8 input parameters: the SO 2 VCD and LH, the 20 surface albedo, the surface height, the O 3 VCD, the solar zenith angle (SZA), the viewing zenith angle (VZA) and the relative azimuth angle (RAA). Table 1 gives a summary of the parameter space. O 3 profiles are classified according to the total column amount of ozone, and the month and latitude zones as specified in the TOMS Version 8 O 3 profile climatology (Bhartia, 2003).
The SO 2 profile plume is taken to have a Gaussian shape, characterized by the total SO 2 VCD loading and centered at a peakconcentration layer height z p , along with a half-width fixed to 2.5 km. In the following, the retrieval of SO 2 layer height refers 25 to the retrieval of the peak-concentration height z p .
Simulations were done on a pressure/temperature/height grid from the US standard atmosphere, with a finer-grid vertical height resolution of 0.25 km below 15 km in order to resolve properly the Gaussian SO 2 plume shape. In total, some 131,072 simulated reflectance spectra have been calculated on a selective parameter grid established by means of a smart sampling technique developed by Loyola et al. (2016). 30 The use of the LIDORT-RRS RT model is necessary, as it enables us to account for the effect of Raman scattering in the atmosphere. Solar irradiances exhibit strong Fraunhofer structures in this part of the UV spectral range, and earthshine spectra are characterized by the "filling-in" of Fraunhofer-solar and telluric-absorber features due to "inelastic" (wavelength-redistributed) rotational-Raman scattering by air molecules. For further details, see Efremenko et al. (2017) and Spurr et al. (2008).
The simulated high-resolution reflectance spectra are then convolved with the TROPOMI Instrument Spectral Response Functions (ISRF) v3.0.0 (released 2018-04-01) 1 . Note that the TROPOMI instrument comprises 450 rows, which are in principle single detectors with their own ISRFs. The signal-to-noise ratio (SNR) is about 1000 in our UV wavelength range. Thus, 5 to account for instrumental noise in the training phase, uncorrelated Gaussian noise with a fixed SNR of 1000 is added to the simulated spectral data.
To extract the information about the layer height and to reduce the dimensionality of the spectral dataset, a principal components analysis (PCA) is applied to the simulated spectra. By thus characterizing the set of simulated measurements with fewer parameters, a simpler, more stable and computationally efficient inversion scheme can be realized. It was found that using 10 10 principal components (PCs) is sufficient to retrieve information about the layer-height. The inclusion of additional PCs beyond 10 did not result in any improvements to LH retrievals, since higher-order PCs are increasingly affected by noise. Figure 6 shows the ratio between the variance of the PCA-derived data set and the total variance of the complete spectra data set (the "explained" variance ratio) as a function of number of PCs included in the PCA. Already with 10 PCs, the error is less than The dimensionality-reduced training spectra, together with the information about the O 3 VCD, the SZA/VZA/RAA angles, and the surface pressure and albedo, are then used as input to train a feed-forward artificial neural network (MultiLayer Perceptron Regression -MLPR), with the SO 2 LH as the output layer. Note here that the SO 2 VCD is not part of the training, since it depends directly on the SO 2 layer height. The training of the MLPR is an iterative process; at each time step, the partial derivatives of the loss function with respect to the model parameters are computed in order to update the parameters. A 20 regularization term is added to the loss function that shrinks model parameters to prevent over-fitting: By building a complex neural network, it is quite easy to perfectly fit the training dataset. When this model is however evaluated on new data (here satellite measurements), it performs very poorly. The regularization modifies the loss function by adding additional terms that penalize large weight vectors and preferring diffuse weight vectors.
After carrying out a PCA and MLPR parameter optimization using closed-loop retrievals to minimize differences between 25 the retrieved and simulated layer heights, the final configuration for the neural network settled on the use of two hidden layers, with 32 nodes in the first layer and 10 nodes in the second.
In the operational phase, the first step is to use the principal component scores acquired during the training phase to transform a given TROPOMI spectral measurement data set to one with a lower dimension. Once this is done, the neural network inverse function is then applied to retrieve the SO 2 LH. 30 In this section, we study the dependency of the layer height retrieval on different parameters. We discriminate between direct (i.e. those parameters affecting the reflectance spectra) and indirect (i.e. affecting the training data and inversion algorithm) dependencies of the retrieved layer height.
-Direct dependencies: Viewing geometry, surface properties, ISRF, Noise level (SNR) 5 -Indirect dependencies: Number of layers in the neural network, number of PCs, parameter ranges for training To proceed, we first divide the training dataset, using 90% of the data for training of the inversion operator and the remaining 10% for testing. Furthermore, we have generated an independent verification dataset modeled with LIDORT_RRS for selected input parameters on a prescribed parameter grid with fixed layer heights of 2.5, 6.5, 13 and 20 km. This verification dataset thus differs from the training dataset, where a "smart" parameter grid was used in place of fixed parameter values.

10
First, we train the FP_ILM operator using 90% of the training dataset (see Table 1), and then apply the trained operator to the remaining 10% training spectra, for which we know the exact SO 2 LH. Figure 1 shows the SO 2 LH difference as a function of solar/viewing geometry, SO 2 VCD, O 3 VCD, albedo and surface pressure. The figure shows clearly a number of marked dependencies in the retrieved layer height, with notably high differences with respect to the real SO 2 LH for low and high layer heights, for low SO 2 VCD well as for high SZA.

15
Regarding the SZA dependency, a cutoff limit of 75 • is set in operational SO 2 retrievals, because at high SZA the light path becomes very long and the noise level increases. Accordingly, we limit the SZA of the spectra used in the training to SZA<70 • in the following.
Clearly, for small SO 2 VCDs, the information content on LH in the spectral signature is very low. It follows that the inclusion of spectra with low SO2 content in the training will have have a negative effect on the entire neural network. In principle for 20 lower SO 2 VCD loadings, more PCs could be included, but at some point also the noise level signature will exceed that of the actual SO 2 absorption.
We have performed several tests in which we limit the training dataset by varying the allowed input parameter ranges. We have finally found an optimal parameter range that allows the retrieval of a broad range of SO 2 LHs even for low SO 2 VCDs. In the training phase, we use only spectra with SO 2 VCD > 20 DU, surface albedo < 0.5 and SZA<70 • to train the final inversion 25 operator for TROPOMI. Figure 2 shows that the error on the retrieved SO 2 LH is less than 2 km and the dependency on SO 2 VCD is reduced.
In addition, we have applied the optimized FP_ILM operator to the independent verification dataset with its prescribed SO 2 LHs. Figure 3 shows the retrieved layer height as a function of SO 2 VCD. As mentioned above, for low SO 2 VCDs, highaltitude layer heights cannot be retrieved -there is always a bias towards low layer heights. Only for SO 2 loadings in excess of 30 20 DU, do we retrieve layer heights with an uncertainty of less than 2 km. To investigate the dependency of the retrieved SO 2 LH as a function of SNR, we have added noise (SNR of 800, 1000 and 1200) to the independent verification spectra. Figure 4 clearly shows that the accuracy of the retrieved layer height improves with increasing SNR.
In the operational retrieval of SO 2 LHs from TROPOMI data, each detector row is effectively a single instrument with its own ISRF. The accuracy of the retrieved SO 2 LH can vary across detector rows. Figure 5 shows the SO 2 LH results when 5 applied to the Sinabung eruption (see Section 4.2). Black dots show the LH for each 50th row, whereas the red cross shows the interpolated LH for the measurement row. Clearly, the retrieved LH is slightly different in each row (within 1 km) due to the different ISRF used for training the NN. Note that in the final operational retrieval, we determine the layer height for a set of fixed detector rows (for which we have trained the FP_ILM separately), and then interpolate the layer height to the actual row.
In this way we avoid jumps in the retrieved layer height between adjacent detector rows. 10 4 Application to TROPOMI data Reflectance spectra from TROPOMI are determined from the L1 solar irradiance and earthshine radiance data (solar irradiance is measured on a daily basis). To correct for Doppler shifts between earthshine and irradiance spectra, we apply the wavelength calibration information from the operational L2 SO 2 product to first calibrate the solar spectrum and then we use the fitted shift and squeeze parameters from the DOAS retrieval to calibrate the earthshine spectra. From this calibrated L1 data, we then 15 calculate reflectances in the wavelength range 310 − 335 nm.
In the following subsections, we have applied the FP_ILM operator to a set of eight major volcanic eruptions measured by TROPOMI. We note that the SO 2 LH retrieval only takes about 2 ms per TROPOMI spectrum, hence even for an extended volcanic plume, the entire LH retrieval can be performed in a matter of seconds, which is important for operational retrieval environments with strong time constraints. To attempt validation of our results, we have performed comparisons with indepen-20 dent MetOp/IASI SO 2 LH data (see Clarisse et al. 2014 for details). The MetOp platforms are on completely different orbits from that of S5P, with widely different overpass times, so that direct satellite comparisons will give only a qualitative validation on the accuracy of retrieved SO 2 LHs from TROPOMI measurements. At the time of writing, IASI data is unfortunately the only source for independent satellite results.

25
A few weeks after the start of the TROPOMI commissioning phase, strong volcanic activity was observed from the Ambae shield volcano (Vanuatu islands). Figure 7 (left) shows the SO 2 VCD retrieved by the standard operational SO 2 algorithm (Theys et al., 2017). We have applied the FP_ILM algorithm to TROPOMI measurements on two occasions (23 November 2017 as well as on 27 July 2018) for overpass times around 02:00h UTC.
On 23 November 2017, two plumes are clearly visible. The first is an aged plume traveling in a south-westerly direction, 30 having a maximum SO 2 VCD of up to 77 DU and a LH of about 10-12 km, while the second plume is moving in a southeasterly direction with SO 2 LH of about 15 km (see Fig. 7). Figure 8 shows the SO 2 LH retrieved for IASI/MetOp-A and -B measurements at overpass times around 11:00h and 21:00h, both events occuring well after the TROPOMI measurements in Fig. 7. Although this timing mis-match makes verification of the SO 2 LH difficult (the plume has already dispersed at the times of the IASI measurements), the maximum SO 2 LH is still about 12-14 km (with a maximum of 20 km), which is in agreement with the FP_ILM results for TROPOMI. The IASI SO 2 LH for the aged plume travelling in the south-west direction is however well below the 10 km LH retrieved by the FP_ILM.

5
On 27 July 2018 a thick SO 2 plume of up to 120 DU was observed, traveling eastwards, see Fig. 9. The retrieved SO 2 LH reached up to 20 km, into the stratosphere, where SO 2 lifetimes are very much longer -here it can form sulfate aerosols and acid rain within weeks. Figure 10 shows the respective IASI SO 2 LH results for overpass times around 13:00h and 23:00h UTC, respectively. The IASI SO 2 LH also peaks at around 20 km, hence showing a very good agreement with the FP_ILM LH results, although the plume has already dispersed.  Figure 13 shows the corresponding SO 2 LH retrieved for IASI/MetOp-A and -B measurements with LH values at about 13 km for the 03:30h UTC overpass and LH up to 18 km around the overpass time of 15:00h UTC; these agree well with the TROPOMI results. Also the Pusat Vulkanologi dan Mitigasi Bencana Geologi 2 (PVMBG, also known as CVGHM) reported at 08:53h UTC a dark gray plume with a high volume of ash that rose at least 20 to 16.8 km. Furthermore the Darwin Volcanic Ash Advisory Center 3 (VAAC) reported that LH values for these Sinabung ash plumes were identified in satellite images, recorded by webcams, and reported by PVMBG continued to rise throughout the day to 13.7 km.

Sierra Negra
On 26 June 2018 a strong eruption at the Sierra Negra shield volcano located on Isabela Island (Galapagos) occurred. This 25 eruption occured in two phases, with an initial very energetic eruption characterized by the opening of five fissures and lava flows, followed by a second event which lasted until 23 August 2018 with decreased gas emissions. TROPOMI was able to measure a very strong SO 2 plume (with loading in excess of 500 DU) only a few minutes after the start of the first eruption. 15), even though the IASI measurements were taken at overpass times of 02:00h and 14:00h UTC (both well before the S5P overpass).

Conclusions
We have developed a new algorithm for the fast and accurate retrieval of SO 2 layer heights from UV earthshine observations of volcanic SO 2 eruptions by the TROPOMI sensor on board the Sentinel-5 Precursor platform. The SO 2 LH retrieval has 5 two phases -(1) a computationally expensive off-line training phase in which the retrieval inverse operator is obtained using the FP_ILM (Full Physics Learning Machine) algorithm; and (2) a very fast operational phase, in which the FP_ILM inverse operator is applied to measured UV reflectance spectra. The FP_ILM combines a Principal Component Analysis with a Neural Network regression using the UV reflectance, O 3 total column, viewing geometries and surface properties as input. Based on an optimized training dataset created with smart sampling techniques, the principal component scores calculated for reflectance 10 spectra in the wavelength range 310-335 nm are used along with the other parameters to train a feed-forward artificial neural network. For S5P/TROPOMI measurement data, an initial dimensionality reduction of the reflectance spectra is performed by applying the PCA-derived Principal Component scores before retrieving SO 2 LH with the trained Neural Network .
The FP_ILM can be used for NRT applications with strict time constraints. S5P/TROPOMI, with its high spatial and spectral resolution, provides a huge amount of data and the computationally intensive direct fitting approaches to SO 2 LH retrieval de-15 veloped so far are not applicable. In contrast, the FP_ILM operator performs SO 2 LH retrieval within about 2 ms per TROPOMI spectrum. Hence even the retrieval on extended volcanic plumes can be performed in a matter of few seconds, allowing for the determination of the SO 2 layer height with an accuracy better than 2 km for SO 2 total column densities larger than 20 DU.
In this paper, we deployed an independent simulated reflectance spectra dataset to investigate the accuracy of SO 2 LH retrievals and their dependencies on a number of different factors and parameters both direct and indirect. In particular, it was 20 found that retrieved SO 2 LH is strongly dependent on the SO 2 VCD for low VCDs <20 DU as well as for high SZAs >75 • .
For high VCDs and low SZA, the SO 2 LH can be retrieved with an accuracy of better than 2 km. We also investigated the dependencies on ISRF and SNR, both of which turned out to be relatively slight effects.
Broad-band spectral scattering and absorption due to sulfate aerosols or volcanic ash plumes will certainly affect SO 2 LH retrievals. Although not considered in the present work, we will address this important issue in a forthcoming paper in this 25 series; aerosols will be accounted for explicitly in the training of the FP_ILM. Nevertheless, we should note that SO 2 and ash are likely to be collocated only for fresh volcanic plumes. For mature plumes, mass differences will ensure that ash and SO 2 plumes are not located at similar altitudes, and the corresponding plumes are thus subject to different wind-direction dispersal.
We have applied the FP_ILM to a number of strong volcanic eruptions observed recently by S5P/TROPOMI. Our SO 2 LH results have been compared to SO 2 LHs retrieved from IR measurements from two MetOp/IASI satellites. Unfortunately the