Lava Volume from Remote Sensing Data: Comparisons with Reverse Petrological Approaches for Two Types of Effusive Eruption

: Five effusive eruptions of Piton de la Fournaise (La R é union) are analyzed to investigate temporal trends of erupted mass and sulfur dioxide (SO 2 ) emissions. Daily SO 2 emissions are acquired from three ultraviolet (UV) satellite instruments (the Ozone Monitoring Instrument (OMI), the Ozone Mapping and Proﬁler Suite (OMPS), and the Tropospheric Monitoring Instrument (TROPOMI)) and an array of ground-based UV spectrometers (Network for Observation of Volcanic and Atmospheric Change (NOVAC)). Time-averaged lava discharge rates (TADRs) are obtained from two automatic satellite-based hot spot detection systems: MIROVA and MODVOLC. Assuming that the lava volumes measured in the ﬁeld are accurate, the MIROVA system gave the best estimation of erupted volume among the methods investigated. We use a reverse petrological method to constrain pre-eruptive magmatic sulfur contents based on observed SO 2 emissions and lava volumes. We also show that a direct petrological approach using SO 2 data might be a viable alternative for TADR estimation during cloudy weather that compromises hot spot detection. In several eruptions we observed a terminal increase in TADR and SO 2 emissions after initial emission of evolved degassed magma. We ascribe this to input of deeper, volatile-rich magma into the plumbing system towards the end of these eruptions. Furthermore, we ﬁnd no evidence of volatile excess in the ﬁve eruptions studied, which were thus mostly fed by shallow degassed magma. to the direct conversion from spectral radiance, and Manual TA the thermodynamic approach.


Introduction
Volcanic eruptions are associated with a multitude of geophysical and geochemical signals, including seismicity (e.g., seismic swarms and volcanic tremor), ground deformation, heat flux, and changes in gas emissions and composition. These signals can be detected and tracked by permanent ground stations and satellite observations [1]. Over the years, development of various technologies for near-real time and continuous data collection has led to improvements in volcano monitoring and tracking of their eruptive behavior. Here, we combine ultraviolet (UV) and infrared (IR) satellite measurements [2][3][4][5][6][7][8] to investigate concurrent trends in lava effusion and sulfur dioxide (SO 2 ) emission during Mapping and Profiler Suite (OMPS) [3,8], and the Tropospheric Monitoring (TROPOMI) [2,7], allow for daily monitoring of volcanic SO2 emissions. Ther (TIR) satellite sensors permit the detection of volcanic hot spots associated activity (e.g., lava flows, lava lakes and domes). Within this category, the Mo lution Imaging Spectroradiometer (MODIS) Level 1B data, provided by NASA's Terra and Aqua satellites, are used as the basis for two automated v itoring systems: Middle InfraRed Observation of Volcanic Activity (MIR MODIS Thermal Alert System (MODVOLC). Further details on the M MODVOLC algorithms are provided by   [6] and Wright Piton de la Fournaise (La Réunion, France) is one of the most active volc world, producing frequent effusive basaltic eruptions of varying duration (av eruption every nine months since 1985 [1,11]), with lava effusion rates typic from one to a few tens of m 3 s -1 [12]. Formed about 4500 years ago, most activity has occurred inside the Enclos Fouqué, a horseshoe-shaped depress eastward to the Indian Ocean. However, eruptions can occur outside the cal ening the surrounding communities. The summit of Piton de la Fournaise composed of two craters: Bory and Dolomieu ( Figure 1). It is also importan despite the intense eruptive activity, during intra-eruptive phases the amou gas is very low [13]. There may be no degassing at all from the volcano betw phases, or the amount can be too low to be measured by satellite instrume amounts are below the detection threshold [13]. This is confirmed by ground urements, which do not detect any SO2 emissions between eruptions [14] demonstrated that the plumbing system at Piton de la Fournaise is compos storage levels connected to each other by sills and dykes ranging from 0.5 to in depth [1,15]. The majority of the eruptions occurring at Piton de la Fourna be fed by the shallowest reservoirs located between 0.5 and 2.5 km beneath crater, with the main reservoir located at about 1.5-2.5 km in depth [15].   [16]. The lava flow outlines and fissure locations are from Chevrel et al. (2021) [17] and OVPF reports ISSN 2610-5101. The high frequency of eruptions at Piton de la Fournaise in recent years provides an extensive dataset to explore the relationships between eruptive heat and gas fluxes. In this study, we estimate lava discharge rates by adapting the methodology of   [18] to five eruptions characterized by different trends in intensity of activity: April 2007, May 2015, August-October 2015, February 2019, and April 2020 (See Figure 1 for the locations of these eruptions). Time-averaged lava discharge rates (TADRs) are used to interpret the evolution of the effusive activity, which in turn is compared with SO 2 emissions detected by the three UV satellite sensors, and the temporal variation of the preeruptive sulfur (S) content in olivine-hosted melt inclusions. The objective of this study is to reconcile estimates of erupted lava volumes from the two MODIS systems (MODVOLC and MIROVA) and the UV sensors used for TADR and SO 2 emission measurements, respectively. The aim is to define the optimum approach for tracking effusive eruptions under different eruptive, meteorological, and instrumental conditions. Furthermore, the analysis of a suite of eruptions at different locations and with different eruptive trends yields insight into the characteristics of the magma reservoir(s) supplying these events [12]. In fact, this comparison allow us to identify different styles of eruptive behavior based on lava and gas emission rates and, in parallel, to confirm that these eruptions were mostly fed by degassed magma resident in the shallow system.

Tropospheric SO 2 Emissions Measured from Space
To quantify eruptive SO 2 emissions, we use daily SO 2 data acquired by the UV OMI, OMPS, and TROPOMI sensors. The SO 2 emissions data are available at the NASA Global Sulfur Dioxide Monitoring website: https://so2.gsfc.nasa.gov/ (last access: 8 October 2021), which collates data from the three instruments. A summary of the main features for each instrument is available in Table A1.

The Ozone Monitoring Instrument
OMI is a hyperspectral UV and visible (VIS) spectrometer capable of detecting and measuring SO 2 associated with volcanic eruptions and passive degassing from space. This instrument is onboard NASA's Aura satellite (operational from September 2004-present), which is in a polar orbit with a local afternoon equator overpass at 13:45. It provides daily, global coverage with a spatial resolution at nadir of 13 × 24 km 2 , but with data gaps since 2009 due to the 'row anomaly' [9]. OMI data products include ozone, SO 2 , and other trace gases such as BrO, HCHO, NO 2 , and OCIO. The UV-2 channel ranges from 307 nm to 383 nm and is used to measure SO 2 with a spectral resolution of 0.45 nm [8,9]. The SO 2 vertical column densities (VCDs) are retrieved using a low-noise principal component analysis (PCA) algorithm [8]. Since the retrieved SO 2 VCD depends on the assumed plume height, the algorithm generates VCDs for three different hypothetical vertical profiles of volcanic SO 2 (3 km, 8 km, and 18 km). If the plume height is known, the correct SO 2 VCD can then be recalculated by linear interpolation. The default SO 2 plume altitude used by NASA's Global Sulfur Dioxide Monitoring website is 8 km (mid-troposphere or TRM), which is appropriate for moderate volcanic eruptions.

The Ozone Mapping and Profiler Suite Instrument
NASA's Suomi NPP/OMPS is a nadir-viewing hyperspectral instrument measuring backscattered UV radiance with a spectral resolution of 1 nm. It is situated in a low Earth orbit, with a local ascending equator overpass at 13:30. OMPS provides daily global coverage with a nadir pixel size of 50 × 50 km 2 . Operational OMPS SO 2 retrievals use the same PCA algorithm used for OMI retrievals [8], and OMPS SO 2 measurements are of similar quality to OMI SO 2 data, but with lower spatial resolution. However, unlike OMI, OMPS currently provides daily global coverage without data gaps. ESA's Sentinel-5 Precursor satellite, also known as Sentinel-5P, carries only one instrument: TROPOMI. This instrument has four spectral channels covering the UV to short-wave infrared (SWIR) wavelength regions. The TROPOMI 310 to 405 nm waveband is used for SO 2 retrievals and has a spectral resolution of 0.54 nm. The Sentinel-5P spacecraft follows a polar orbit with a local equator crossing at 13:30 (ascending node). With a spatial resolution of 7 × 3.5 km 2 , TROPOMI provides daily SO 2 measurements with higher spatial resolution than OMI and OMPS [2]. TROPOMI SO 2 retrievals employ the differential optical absorption spectroscopy (DOAS) method using a combination of three fitting windows, namely 312-326 nm as standard, plus 325-335 nm or 360-390 nm as alternatives to avoid signal saturation for high SO 2 loading [7].

NOVAC Network
As well as SO 2 data from satellite instruments, we use ground-based measurements from a set of three DOAS instruments from the Network for Observation of Volcanic and Atmospheric Change (NOVAC) operated by the Observatoire Volcanologique du Piton de la Fournaise (OVPF) of the Institut de Physique du Globe de Paris. The stations have been scanning over the Dolomieu crater's rim since September 2007, i.e., after the April 2007 eruption. The stations are at distances of 4 to 6 km W of the crater (Figure 1), and cover plumes transported north, west, or south from the volcano. NOVAC measurements involve the acquisition of UV spectra (280-420 nm) of skylight over a conical scanning plane enclosing the volcano. The spectra obtained are analyzed using the DOAS technique [19] to derive SO 2 columns in the range of 310-325 nm. Plume height and direction are obtained from triangulation of measurements from two stations, while plume speed is obtained from meteorological stations or atmospheric transport models. For this study, we used wind speed from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 re-analysis database (0.25 × 0.25 deg, summit level, hourly). The typical time resolution of a gas flux measurement is 5 min, and the uncertainty of a fully scanned plume lies at about 30% [19]. However, the possibility and quality of the measurements depend on the wind direction, the weather conditions, and the magnitude of the eruption, as well as the vent location. Gas emissions escape detection by NOVAC instruments during nighttime, when transported to the east of the volcano, when the plume is not lifted above the level of the Dolomieu crater's rim, or when its altitude is several (>5) km high above the station. Incompletely scanned plumes would lead to underestimation of the flux, but the spatial information could be used to determine the direction and altitude of the plume.

Materials and Methods
We processed OMI, OMPS, and TROPOMI SO 2 data for all eruptions at Piton de la Fournaise since 2004, deriving and comparing daily SO 2 emissions. Figure A1 gives the time-series of SO 2 mass during the selected eruptions. In parallel, TADR data were acquired for the MODVOLC and MIROVA systems and estimated from manual processing of the MODIS images. By integrating these data through an eruption, the cumulative lava volume can be determined and compared to field observations, assuming the latter are the most accurate, to validate the satellite-based estimates. In addition, total SO 2 emissions for each eruption are estimated for a range of sulfur (S) contents within melt inclusions and matrix glasses using a direct petrological approach. A reverse petrological approach is also used to estimate the daily pre-eruptive magmatic sulfur content during the eruptions by fixing the residual sulfur content within the lava matrix and assuming that the observed SO 2 emissions are derived from the satellite-based lava mass.

Field Measurements
Generally, field measurements only include the bulk lava volumes erupted. However the total lava volume should also include the cone and pyroclasts emitted during the event. A linear correlation was established by Michon et al. (2013) [20] between the lava Remote Sens. 2022, 14, 323 5 of 26 and pyroclast volume, allowing for the estimation of this missing volume. Usually it represents about 10% of the total volume or less. The field bulk lava volume data presented in this study came from the analysis of Derrien (2019) [21]. The application of digital image processing techniques to (aerial) photographs (a.k.a. photogrammetry) associated with 3D models provide a fine analysis of lava flows volumes that takes into account the morphological evolution of the ground. The reader may refer to Derrien (2019) [21] for further information on the techniques used to acquire this data. In the following sections, we also discuss the Dense Rock Equivalent (DRE) lava volume, which is obtained by removing the vesicularity and crystallinity of the lava from the bulk lava volume.

MODIS Image Processing
TADRs for Piton de la Fournaise eruptions were calculated from thermal data acquired by two automated systems: MIROVA and MODVOLC. For the MODVOLC approach, the data are publicly available on http://modis.higp.hawaii.edu/ (last access: 8 October 2021). MIROVA data were provided by the University of Turin. Both MODVOLC and MIROVA use MODIS Level 1b data collected from NASA's Terra and Aqua satellites [4][5][6]. However, MODVOLC provides thermal anomaly data in terms of total lava flow radiance (L f low ), whereas MIROVA gives those data in terms of Volcanic Radiative Power (VRP = [1.89 × 10 7 ] × L f low ) [22,23]. Terra and Aqua MODIS acquire data in 36 spectral bands ranging from 0.4 to 14.4 µm and provide global coverage every one to two days following the EOS orbit. The bands of interest have a nominal resolution of 1 km at nadir [5]. Our goal was to compare the TADRs derived from MODVOLC and MIROVA and to validate the TADRs using manual analysis of MODIS L1b data (hereafter referred to as "Manual"). Note that the TADR and the effusive trend are provided by considering only the MODIS Level 1b data acquired during good weather conditions. Hence the raw L1b MODIS radiance data were downloaded in order to manually select the hot spots and background pixels using the ENVI software. As part of this manual approach, spectral radiances were extracted from four spectral bands of interest:
Following pixel selection, the pixel spectral radiances were corrected for surface emissivity, atmospheric transmission, atmospheric emission, and surface reflection effects using the following equation: where L is the Planck function, L(λ, T S ) is the surface radiance for a surface at temperature T S measured at wavelength λ, L(λ, T * ) is the at-sensor radiance, ε λ is spectral emissivity, τ λ is the atmospheric spectral transmissivity, and L u (λ) is atmospheric upwelling radiance. The three last parameters were obtained using MODTRAN software [22]. Further details of the procedure are given in Harris (2013) [24]. The most important factor is saturation, which happens when the amount of emitted radiance that the MODIS sensor can detect is exceeded. To determine the pixel fraction occupied by a hot target necessary to saturate the different bands, Equation (2) is used [25]: where L(T b ) is the background radiance for the cold surface at temperature T b , and L(T hot ) is the radiance of the hot spot. For a given hot spot and background temperature, we were able to determine p, the pixel fraction needed to reach the saturation level of the sensor, L(T sat ) with T sat being the saturation temperature. The saturation temperatures were fixed at 60 • C, 180 • C, and 130 • C for Bands 22, 21, and 31 and 32, respectively [25]. Once the saturation point is reached, only a single value is recorded, so even if the area and temperature of the hot spot increases, the brightness temperature will remain the same. Thus, any derived value (e.g., TADR) using saturated values will provide a minimum limit on the value, where Figure 2 shows the upper limit of the size/temperature of the feature that can be measured before this occurs for a T b of 25 • C. By taking a 1 km 2 pixel (e.g., at nadir) and a hot spot temperature of 500 • C, Figure 2 illustrates that only a small portion of a pixel (3.50%) is needed to complete a saturation in the MIR compared to those in TIR (12.35-13.30%). It is also important to point out that Band 22 reaches the saturation level for a much smaller area (1200 m 2 ) than Band 21 (35,000 m 2 ).
sensor, ( ) with being the saturation temperature. The saturation were fixed at 60 °C, 180 °C, and 130 °C for Bands 22, 21, and 31 and 32, res Once the saturation point is reached, only a single value is recorded, so ev and temperature of the hot spot increases, the brightness temperature w same. Thus, any derived value (e.g., TADR) using saturated values will pr mum limit on the value, where Figure 2 shows the upper limit of the size/t the feature that can be measured before this occurs for a of 25 °C. By t pixel (e.g., at nadir) and a hot spot temperature of 500 °C, Figure 2 illustra small portion of a pixel (3.50%) is needed to complete a saturation in the M to those in TIR (12.35-13.30%). It is also important to point out that Band saturation level for a much smaller area (1200 m 2 ) than Band 21 (35,000 m 2 )

Estimation of Lava Discharge Rate from MODIS Data
Two methods were applied in this study to estimate TADR that use a b lava volume. The first method is a linear relationship directly relating the sa spectral radiance and TADR (Direct Conversion, DC): where is the hot spot radiance minus the background radiance su relevant MODIS pixels, and is a constant equal to 0.128 m 5 J -1 , as determin May-July 2003 eruption at Piton de la Fournaise [26,27]. The second way TADR is to utilize a thermodynamic approach (TA) [28]: where is the Stefan Boltzmann constant (5.67 × 10 W m -2 K -4 ), ℎ is t heat transfer coefficient (10 W m -2 K -1 , [24]), and are the temperature surface and background, respectively, is the bulk density of the lava ( the bulk specific heat (J kg -1 K -1 ), is the latent heat of crystallization (J k cooling range between the vent and the flow front [26] (K), is the crystalli ing through ∆ , and is the active lava area (m 2 ). Equation (4) can be linear relationship, where is an empirical parameter that converts = ×

Estimation of Lava Discharge Rate from MODIS Data
Two methods were applied in this study to estimate TADR that use a bulk (not DRE) lava volume. The first method is a linear relationship directly relating the satellite-derived spectral radiance and TADR (Direct Conversion, DC): where L λ− f low is the hot spot radiance minus the background radiance summed for all relevant MODIS pixels, and c is a constant equal to 0.128 m 5 J -1 , as determined during the May-July 2003 eruption at Piton de la Fournaise [26,27]. The second way to obtain the TADR is to utilize a thermodynamic approach (TA) [28]: where σ is the Stefan Boltzmann constant (5.67 × 10 −8 W m -2 K -4 ), h c is the convective heat transfer coefficient (10 W m -2 K -1 , [24]), T c and T a are the temperature (K) of the lava surface and background, respectively, ρ m is the bulk density of the lava (kg m -3 ), c p is the bulk specific heat (J kg -1 K -1 ), c L is the latent heat of crystallization (J kg -1 ), ∆T is the cooling range between the vent and the flow front [26] (K), ϕ is the crystallization in cooling through ∆T, and A lava is the active lava area (m 2 ). Equation (4) can be reduced to a linear relationship, where m b is an empirical parameter that converts A lava to TADR. The constants m and b combine all the assumed values that are the same for all measurements [27,29], where m = 2.07 × 10 4 and b = 6.12 × 10 7 . Using a value of 3.38 × 10 −4 for m b , Equation (5) provides nearly the same results as the direct conversion from spectral radiance (Equation (3)), except for the April 2007 eruption, where a m b value of 3.53 × 10 −4 gives the best fit. In addition, following Coppola et al. (2009) [27], for the April 2007 eruption, the maximum TADR is increased by 30% due to the flux lost to the ocean. By integrating TADR through the whole duration of the eruption, we acquire the volume of lava erupted with 30% uncertainty [22,29].
Instead, MIROVA data are provided in terms of Volcanic Radiative Power (VRP in W). Hence, TADR is calculated using a coefficient called radiant density (c rad in J m -3 ), which is based on the silica content of the erupted lavas, as follows [22,30]:

Petrologic Method
The petrological approach assumes that melt inclusions trapped within olivine phenocrysts in the magma chamber represent the pre-eruptive volatile contents dissolved in the magma and that no accumulation of exsolved volatiles occurs in the reservoir. We assume that there are no additional contributions from independent fluid and solid phases, and that melt inclusions record variable degrees of magma degassing in the plumbing system [31]. As a consequence, the source of the degassed S is only the melt phase of the magma. The direct (standard) petrologic method, therefore, estimates the mass of SO 2 (M SO 2 ) released as Results are given with 45 % uncertainty. V lava is the total (bulk) volume of lava emitted during the eruption (m 3 -here the MIROVA data were used), c MI S and c matrix S are the sulfur concentrations recorded in the melt inclusion and the degassed matrix (ppm), respectively, ρ m is the average bulk melt density assuming a melt density of 2940 kg m -3 , φ is vesicularity, ε crx is the average volume fraction of phenocrysts in the magma [32], and α (kg) is a constant clustering all the parameters cited above. Table 1 summarizes the crystal fraction, vesicularity, and the bulk density data used for the studied eruptions. Table 1. Summary of petrologic method data used to estimate the sulfur content and SO 2 mass, where ε crx is the crystalinity and φ the vesicularity (http://wwwobs.univ-bpclermont.fr/SO/televolc/ dynvolc/ (last access: 14 June 2021)). According to petrological analysis [33], c matrix S can be fixed at 160 ± 60 ppm during high intensity phases (e.g., lava fountains). This residual value can increase up to 230 ± 30 ppm during the phases of fast magma ascent and incomplete degassing [13,33] parameter and the daily SO 2 and lava masses detected from space, Equation (7) may be arranged to determine the pre-eruptive magma S content as Vesicularity and crystallinity, included in the α parameter, can vary during eruptions (see Table 1) and so are adjusted accordingly. In this way, a time series of the pre-eruptive sulfur content can be estimated with 50 % uncertainty. In addition, we assume that the SO 2 is not significantly dissolved in condensed water or absorbed on ash and/or Pele's hairs in the volcanic plume.

Comparisons with NOVAC
SO 2 masses calculated from the satellite data were first compared with the SO 2 flux measured from ground-based NOVAC stations for the five selected eruptions ( Figure 3). Despite some differences, Figure 3 shows good agreement for the eruptions occurring in 2015, where the same trends were apparent in all datasets (OMPS, OMI, and NOVAC). However, for the February 2019 and April 2020 eruptions, the NOVAC fluxes were much lower than the satellite measurements. Vesicularity and crystallinity, included in the parameter, can vary during eruptions (see Table 1) and so are adjusted accordingly. In this way, a time series of the preeruptive sulfur content can be estimated with 50 % uncertainty. In addition, we assume that the SO2 is not significantly dissolved in condensed water or absorbed on ash and/or Pele's hairs in the volcanic plume.

Comparisons with NOVAC
SO2 masses calculated from the satellite data were first compared with the SO2 flux measured from ground-based NOVAC stations for the five selected eruptions ( Figure 3). Despite some differences, Figure 3 shows good agreement for the eruptions occurring in 2015, where the same trends were apparent in all datasets (OMPS, OMI, and NOVAC). However, for the February 2019 and April 2020 eruptions, the NOVAC fluxes were much lower than the satellite measurements.

MODIS Analysis
A comparison can be made between the lava volume obtained by the integrated MODIS-derived TADR with the volume measured in the field. For both discharge rates and volumes, the estimations correspond to those derived from the MIROVA and MODVOLC automatic systems and the Manual method. Note that the MIROVA and MODVOLC systems do not consider the atmospheric correction, contrary to the Manual approaches. To simplify, in the following sections, we only consider the first Manual approach (Equation (3)). All comparisons are shown in Figure 4.

MODIS Analysis
A comparison can be made between the lava volume obtained by the integrated MODIS-derived TADR with the volume measured in the field. For both discharge rates and volumes, the estimations correspond to those derived from the MIROVA and MODVOLC automatic systems and the Manual method. Note that the MIROVA and MODVOLC systems do not consider the atmospheric correction, contrary to the Manual approaches. To simplify, in the following sections, we only consider the first Manual approach (Equation (3)). All comparisons are shown in Figure 4.

April 2007
The April 2007 eruption produced a total of 185 × 10 6 m 3 of lava (bulk volume) [21]. As the eruption was very intense, instruments were saturated, with there being between 3 and 10 saturated pixels per image (i.e., 17-60% of the anomalous pixels) so that minimum bounds on TADR were all that could be given. Figure 4a shows a peak on 5 April that was recorded by all methods. On the MIROVA trend, a peak of up to 200 m 3 s -1 is obtained, which differs significantly from peak values obtained with the other methods (Manual: 48 m 3 s -1 , MODVOLC: 51 m 3 s -1 ). Note that this value (200 m 3 s -1 ) was manually fixed, and so was not derived from MIROVA data [27]. This peak occurred a few hours before the beginning of the collapse of the Dolomieu crater, which occurred during the night of 5 April and lasted several days. We note that the lava discharge rate was nearly Note that on the MIROVA dataset, the peak on 6 April 2007 that was manually corrected (200 m 3 s -1 ) was taken into account. By removing this corrected value, the lava volume estimation decreases to 52.1 ± 26.0 × 10 6 m 3 . The uncertainty is 30% except for the April 2007 eruption, where it rises to 50% due to saturation issues.

April 2007
The April 2007 eruption produced a total of 185 × 10 6 m 3 of lava (bulk volume) [21]. As the eruption was very intense, instruments were saturated, with there being between 3 and 10 saturated pixels per image (i.e., 17-60% of the anomalous pixels) so that minimum bounds on TADR were all that could be given. Figure 4a shows a peak on 5 April that was recorded by all methods. On the MIROVA trend, a peak of up to 200 m 3 s -1 is obtained, which differs significantly from peak values obtained with the other methods (Manual: 48 m 3 s -1 , MODVOLC: 51 m 3 s -1 ). Note that this value (200 m 3 s -1 ) was manually fixed, and so was not derived from MIROVA data [27]. This peak occurred a few hours before the beginning of the collapse of the Dolomieu crater, which occurred during the night of 5 April and lasted several days. We note that the lava discharge rate was nearly zero on 6 April. A second peak was recorded after the next and smaller collapse of the Dolomieu crater occurring on 12 April. TADR estimation from the MODVOLC system (78 m 3 s -1 ), MIROVA (70 m 3 s -1 ), and the Manual method (73 m 3 s -1 ) were in good agreement. The next few days were marked by a decline in the discharge rates, occasionally interrupted by minor fluctuations at the end of April (Figure 4a). Assuming that the bulk field volume is accurate, the manual and MODVOLC approaches largely underestimate the volume (Table 2). However, MIROVA gave a higher estimate closer to that expected by using the corrected dataset, i.e., using the peak in TADR on 6 April of 200 m 3 s -1 . Without considering this correction, MIROVA data provided a lava volume of 52.1 ± 26.0 × 10 6 m 3 , which matched with MODVOLC and the Manual results. In Figure 5, one may observe that the time-series of seismically-derived TADR from Duputel and Rivera (2019) [34] fit with the MIROVA data despite their higher values (up to 600 m 3 s -1 ) [34]. zero on 6 April. A second peak was recorded after the next and smaller collapse of the Dolomieu crater occurring on 12 April. TADR estimation from the MODVOLC system (78 m 3 s -1 ), MIROVA (70 m 3 s -1 ), and the Manual method (73 m 3 s -1 ) were in good agreement. The next few days were marked by a decline in the discharge rates, occasionally interrupted by minor fluctuations at the end of April (Figure 4a). Assuming that the bulk field volume is accurate, the manual and MODVOLC approaches largely underestimate the volume (Table 2). However, MIROVA gave a higher estimate closer to that expected by using the corrected dataset, i.e., using the peak in TADR on 6 April of 200 m 3 s -1 . Without considering this correction, MIROVA data provided a lava volume of 52.1 ± 26.0 ×10 6 m 3 , which matched with MODVOLC and the Manual results. In Figure 5, one may observe that the time-series of seismically-derived TADR from Duputel and Rivera (2019) [34] fit with the MIROVA data despite their higher values (up to 600 m 3 s -1 ) [34].

May 2015
This event was characterized by a TADR peak of 30 m 3 s -1 at the beginning of the eruption on May 17-18 followed by an exponential decrease (Figure 4a). Despite good agreement between the approaches, we observe a significant discrepancy on May 17 (MODVOLC: 29 m 3 s -1 , MIROVA: 25 m 3 s -1 , Manual: 14 m 3 s -1 ). In terms of cumulative volume, the high lava discharge rates from MODVOLC led to the highest cumulative volume (8.9 × 10 6 m 3 ). These were greater than MIROVA (6.0 × 10 6 m 3 ), but MIROVA matched the lava volume measured in the field (5.7 × 10 6 m 3 ) (Table 2, Figure 4b). Figure 4a shows a short-lived TADR peak of 50 m 3 s -1 at the beginning of the eruption. The activity started on August 24 with high TADRs along a fissure system, but some differences resulted from the distinct methods (MIROVA: 50 m 3 s -1 , Manual: 24 m 3 s -1 ). Following this initial peak, TADRs remained nearly constant (5-10 m 3 s -1 up to 13-20 m 3 s -1 for MODVOLC data) during the second phase of the eruption that lasted from September 11 until October 16. The last phase involved three final pulses each interrupted by 4 days of inactivity (Figure 4a) Table 2). However, the MODVOLC algorithm overestimated the lava volume emitted (73.3 ×10 6 m 3 ) during the August-October eruption, an event that was particularly long-lived, as it lasted 65 days, resulting in a cumulative over-estimation.

February 2019
For this eruption, TADR was steady over time and below 10 m 3 s -1 , except at the end when TADR increased to 37 m 3 s -1 according to MIROVA (Figure 4a). We identified a peak at the end of the eruption that was not as high in the manually processed data as in the MIROVA data (19 m 3 s -1 versus 37 m 3 s -1 ; Figure 4a). In general, MODVOLC-derived values were again much higher than all other methods, ranging from 9 to 31 m 3 s -1 . In addition, the MODVOLC system stopped detecting hot spots after March 2. Hence, no estimation of TADR could be made for the end of February, thus missing the intense terminal phase of the eruption. MIROVA and field estimates of erupted lava volumes were similar, at 13.6 × 10 6 m 3 and 14.5 × 10 6 m 3 , respectively ( Table 2). The MODVOLC and manual approaches provided cumulative volumes that were also similar, but lower, at 9.8 and 9.4 ×10 6 m 3 , respectively ( Table 2).

April 2020
Only three detections were made by MODVOLC. This was due to the extremely cloudy conditions during this eruption, leading to poor detection rates. We saw two TADR peaks on 4 and 5 April (Figure 4a). Although MIROVA and manual values were generally similar (with a difference of 6 m 3 s -1 ), MIROVA estimated a much higher TADR on 5 April (22 against 8 m 3 s -1 for Manual). In terms of cumulative volume, we obtained a volume that was between 2 and 5 × 10 6 m 3 less than the bulk field volume ( Table 2).

Summary of Thermally-Derived TADR Data
Although consistent in terms of trend, there was some divergence in the absolute values derived from each approach (Figure 4). Discrepancies arise from variable detection limits of the algorithms employed by the automatic systems, as well as saturation problems (Table 3), plus human error regarding the manual method. In addition, the MIROVA and MODVOLC systems do not correct for emissivity, atmospheric transmission, and surface reflection. Hence, the actual surface-leaving spectral radiances were underestimated by 25% due to atmospheric absorption. However, despite this, in general, MODVOLC seemed to overestimate TADR. Instead, we obtained a good correlation between MIROVA and Manual data ( Figure 6). Table 3. Summary of the total pixels for each eruption counted by each method. Note that the MIROVA algorithm resampled MODIS images into 1 km equal area pixels (e.g., high zenith angle). Hence, for some cases, one original MODIS pixel may be divided into several "MODIS-MIROVA" pixels. That explains the differences in total pixels between MIROVA and the other approaches. "#/% pixels sat" = number/percentage of saturated pixel.  Table 3. Summary of the total pixels for each eruption counted by each method. Note that the MI-ROVA algorithm resampled MODIS images into 1 km equal area pixels (e.g., high zenith angle). Hence, for some cases, one original MODIS pixel may be divided into several "MODIS-MIROVA" pixels. That explains the differences in total pixels between MIROVA and the other approaches. "#/% pixels sat" = number/percentage of saturated pixel.  By integrating TADR through time, all methods provided similar bulk cumulative lava volumes. However, assuming that the bulk field volumes are accurate, we find that in April 2020 and April 2007, all approaches significantly underestimated the erupted lava volumes ( Table 2). This was due to bad weather and/or saturation of the sensor, which caused data loss and under-estimation. However, the total lava volumes for the year 2015 and February 2019 were in good agreement with the ground truth, except for those given by the MODVOLC system, which overestimated despite missing days (Table 3). We address the cause for this in the discussion. Figure 7 shows the retrieved magmatic sulfur contents using the reverse petrological approach and the erupted lava masses obtained from MODIS TIR satellite data. By integrating TADR through time, all methods provided similar bulk cumulative lava volumes. However, assuming that the bulk field volumes are accurate, we find that in April 2020 and April 2007, all approaches significantly underestimated the erupted lava volumes ( Table 2). This was due to bad weather and/or saturation of the sensor, which caused data loss and under-estimation. However, the total lava volumes for the year 2015 and February 2019 were in good agreement with the ground truth, except for those given by the MODVOLC system, which overestimated despite missing days (Table 3). We address the cause for this in the discussion. For the April 2007 eruption, all SO2 masses estimated using the direct petrological method were lower than the 270 kt of SO2 detected by the OMI sensor, indicating a possible sulfur "excess" during this eruption [36,37]. The average pre-eruptive S content in melt inclusions obtained with the reverse petrological method by fixing at 160 ppm during high flux phases or 230 ppm during low degassing (corresponding to phases of intense fountains) was very high (up to 30,000 ppm). We expected much lower values (ca. 1050 ppm; [36]) according to the measured S within melt inclusions from petrological analysis ( Figure 8). However, by using the time-series of seismically-derived TADR from Duputel and Rivera (2019) [34], the total SO2 mass estimated from the direct petrological approach was about 370 kt. Using the lava masses obtained from their data and the reverse petrological approach, we estimated an average S within melt inclusion between 1050 and 1250 ppm (Figure 7a), fully consistent with petrological analysis [36]. (a-e) Estimation of sulfur content within melt inclusions using the reverse petrological approach. The larger the symbol, the better the match with the total SO 2 mass detected by the UV satellite sensors. According to petrological analysis, the average residual sulfur content in the matrix can be fixed at 160 ppm [33]. This is represented by the black line. Note that for the April 2007 eruption, the sulfur content estimation used the seismically-derived TADR data from Duputel and Rivera (2019) [34] and the three points corresponding to the following S couples (MI-Matrix: 1050-160 ppm; 1150-260 ppm; 1250-360 ppm), provide exactly the total SO 2 mass estimated from the direct petrological approach. Results are given with 45% uncertainty. For the April 2007 eruption, all SO 2 masses estimated using the direct petrological method were lower than the 270 kt of SO 2 detected by the OMI sensor, indicating a possible sulfur "excess" during this eruption [36,37]. The average pre-eruptive S content in melt inclusions obtained with the reverse petrological method by fixing c matrix S at 160 ppm during high flux phases or 230 ppm during low degassing (corresponding to phases of intense fountains) was very high (up to 30,000 ppm). We expected much lower values (ca. 1050 ppm; [36]) according to the measured S within melt inclusions from petrological analysis ( Figure 8). However, by using the time-series of seismically-derived TADR from Duputel and Rivera (2019) [34], the total SO 2 mass estimated from the direct petrological approach was about 370 kt. Using the lava masses obtained from their data and the reverse petrological approach, we estimated an average S within melt inclusion between 1050 and 1250 ppm (Figure 7a), fully consistent with petrological analysis [36]. emote Sens. 2022, 13, x FOR PEER REVIEW Figure 8. Petrological analysis of sulfur in melt inclusions at Piton de la Fournaise fro [13,38]. Inclusions having a sulfur content of ca. 1050 ppm are considered to come fro reservoir. Those that have a higher sulfur content (>1150 ppm) are related to deepe while partly degassed and shallow seated magmas are recorded by S contents <1000 Given low pre-eruptive S contents of between 100 and 400 ppm, estima sions for the May 2015 eruption were consistent with an eruption largely fed magma. Both OMI and OMPS SO2 data imply ~200-500 ppm of sulfur in m (Figure 7b). Similar values were obtained for the August-October 2015 e OMPS data. OMI data yielded slightly higher values of ~300-650 ppm S (F despite this slight difference, the estimated sulfur contents were in good a tween the two sensors.

Magmatic Sulfur Content Estimations
However, for the February 2019 eruption, the difference between the sensors was greater. Whereas TROPOMI and OMI indicated a similar rang sulfur content (300-800 ppm), OMPS gave a higher range (800 to 1600 ppm that deep, undegassed magma [12,27,33]was also involved in this eruptio The situation was similar for the April 2020 eruption, where TROPOMI around 2000 ppm S and OMI between 900 and 1600 ppm S (Figure 7e). Figu show the OMPS results, as no magmatic sulfur content values matched the SO2 detected by the sensor. On 4 and 6 April, TROPOMI and OMPS meas 26.4 kt of SO2 ( Figure A1), respectively, whereas OMI data gaps led to an un of the total SO2 emissions during the eruption. Figure 9 shows the daily pre-eruptive magmatic sulfur contents estimate of the reverse petrological approach using satellite-derived SO2 masses and undegassed sulfur content in the eruptive products. The April 2007 eruptio sented, as the values were greater than 30,000 ppm S (using the MIROVA which is unrealistic. Note that this approach is dependent on the detection of sensors. Consequently, if a sensor detected no SO2, we underestimated the em added an error to the estimated pre-eruptive sulfur content. Nevertheless, obtained time-series that were in good agreement with each other (Figure 9).  [13,38]. Inclusions having a sulfur content of ca. 1050 ppm are considered to come from the sea level reservoir. Those that have a higher sulfur content (>1150 ppm) are related to deeper mafic inputs, while partly degassed and shallow seated magmas are recorded by S contents <1000 ppm.

Temporal Variation of Sulfur Content
Given low pre-eruptive S contents of between 100 and 400 ppm, estimated SO 2 emissions for the May 2015 eruption were consistent with an eruption largely fed by degassed magma. Both OMI and OMPS SO 2 data imply~200-500 ppm of sulfur in melt inclusions (Figure 7b). Similar values were obtained for the August-October 2015 eruption using OMPS data. OMI data yielded slightly higher values of~300-650 ppm S (Figure 7c), but despite this slight difference, the estimated sulfur contents were in good agreement between the two sensors.
However, for the February 2019 eruption, the difference between the three satellite sensors was greater. Whereas TROPOMI and OMI indicated a similar range of magmatic sulfur content (300-800 ppm), OMPS gave a higher range (800 to 1600 ppm), suggesting that deep, undegassed magma [12,27,33] was also involved in this eruption (Figure 7d). The situation was similar for the April 2020 eruption, where TROPOMI data indicates around 2000 ppm S and OMI between 900 and 1600 ppm S (Figure 7e). Figure 7e does not show the OMPS results, as no magmatic sulfur content values matched the total mass of SO 2 detected by the sensor. On 4 and 6 April, TROPOMI and OMPS measured 12.1 and 26.4 kt of SO 2 ( Figure A1), respectively, whereas OMI data gaps led to an underestimation of the total SO 2 emissions during the eruption. Figure 9 shows the daily pre-eruptive magmatic sulfur contents estimated on the basis of the reverse petrological approach using satellite-derived SO 2 masses and by fixing the undegassed sulfur content in the eruptive products. The April 2007 eruption is not represented, as the values were greater than 30,000 ppm S (using the MIROVA lava volume), which is unrealistic. Note that this approach is dependent on the detection of SO 2 by the UV sensors. Consequently, if a sensor detected no SO 2 , we underestimated the emitted SO 2 and added an error to the estimated pre-eruptive sulfur content. Nevertheless, in general, we obtained time-series that were in good agreement with each other (Figure 9). Remote Sens. 2022, 13, x FOR PEER REVIEW 15 of 26 Results were still lower than the value of 1050 ppm S expected for the sea-level reservoir, which was in good agreement with the general view of involvement of variably degassed magmas in many recent eruptions [33] (Figure 9). For the February 2019 eruption, was, on average, equal to 521 ± 285 ppm, which was obtained using OMI and TROPOMI data. However, OMPS data indicated a melt enriched in sulfur, which could be explained by the visible peak at the end of the eruption (Figure 9c). A similar observation could be made for the August-October 2015 eruption (Figure 9e). Despite an almost constant value of 283 ± 191 ppm (according to OMI data), some peaks were observed, notably during the third and last phases of the eruption (16 October-2 November), suggesting the involvement of a deep undegassed melt.

Daily Volume Estimation from Sulfur Content
We also re-arranged Equation (8) to calculate the expected DRE lava volume using the measured daily SO2 emissions and fixing (160 ppm) and . The first approach was to fix at 1050 ppm corresponding to the sea level magma reservoir. The second method used the above (Figure 9). Values higher than 2000 ppm were replaced by 1250 or 1600 ppm, as observed in the petrological analysis ( Figure 8). In this case, the DRE field lava volume was used as a reference.

Fixed Sulfur Content within Melt Inclusions
By fixing at 1050 ppm, we generally obtained lower lava volumes than those derived from the field data for the two eruptions in 2015 (May 2015, Aug-Oct 2015) ( Table  4). OMI-derived total lava volumes were, in general, lower relative to field data (Table 4). OMPS data overestimated the lava volume for the April 2020 and February 2019 eruptions but underestimated lava volumes for the eruptions in 2015. TROPOMI overestimated the field volume, being 8.4 × 10 m 3 and 3.3 × 10 m 3 , respectively, (Table 4) for the April 2020 eruption. In addition, in February 2019 the field volume was slightly higher than our estimates (9.4 × 10 m 3 versus 6.3 × 10 m 3 ) (Table 4). Consequently, fixing the pre- Results were still lower than the value of 1050 ppm S expected for the sea-level reservoir, which was in good agreement with the general view of involvement of variably degassed magmas in many recent eruptions [33] (Figure 9). For the February 2019 eruption, c MI S was, on average, equal to 521 ± 285 ppm, which was obtained using OMI and TROPOMI data. However, OMPS data indicated a melt enriched in sulfur, which could be explained by the visible peak at the end of the eruption (Figure 9c). A similar observation could be made for the August-October 2015 eruption (Figure 9e). Despite an almost constant c MI S value of 283 ± 191 ppm (according to OMI data), some peaks were observed, notably during the third and last phases of the eruption (16 October-2 November), suggesting the involvement of a deep undegassed melt.

Daily Volume Estimation from Sulfur Content
We also re-arranged Equation (8) to calculate the expected DRE lava volume using the measured daily SO 2 emissions and fixing c matrix S (160 ppm) and c MI S . The first approach was to fix c MI S at 1050 ppm corresponding to the sea level magma reservoir. The second method used the c MI S above (Figure 9). Values higher than 2000 ppm were replaced by 1250 or 1600 ppm, as observed in the petrological analysis ( Figure 8). In this case, the DRE field lava volume was used as a reference.

Fixed Sulfur Content within Melt Inclusions
By fixing c MI S at 1050 ppm, we generally obtained lower lava volumes than those derived from the field data for the two eruptions in 2015 (May 2015, Aug-Oct 2015) ( Table 4). OMI-derived total lava volumes were, in general, lower relative to field data ( Table 4). OMPS data overestimated the lava volume for the April 2020 and February 2019 eruptions but underestimated lava volumes for the eruptions in 2015. TROPOMI overestimated the field volume, being 8.4 × 10 6 m 3 and 3.3 × 10 6 m 3 , respectively, (Table 4) for the April 2020 eruption. In addition, in February 2019 the field volume was slightly higher than our estimates (9.4 × 10 6 m 3 versus 6.3 × 10 6 m 3 ) (Table 4). Consequently, fixing the pre-eruptive sulfur content at 1050 ppm may be unreasonable, and a lower value may be more appropriate (see next section). Table 4. Cumulative DRE volumes of erupted lava (×10 6 m 3 ) estimated from field measurements, MIROVA, and the reverse petrological approach using OMI, OMPS, and TROPOMI SO 2 data. "OMI/OMPS/TROP 1050" = pre-eruptive sulfur content within melt inclusions fixed at 1050 ppm; "OMI/OMPS/TROP t" = pre-eruptive sulfur content varies over time. Note that the MIROVA cumulative lava volume is expressed as a bulk volume and not a DRE volume.

Eruption
Cumulative

Variable Pre-Eruptive Sulfur Contents
As no realistic predictions for the sulfur content for the April 2007 eruption were obtained using the MIROVA data, this method could not be applied for this event.
Results showed similar daily lava volumes for April 2020 compared to the previous method (Figure 10a,b). For the February 2019 eruption, the cumulative volume derived using TROPOMI (11.6 × 10 6 m 3 ) and OMI (9.2 × 10 6 m 3 ) were in good agreement with field data (9.4 × 10 6 m 3 ). For the August-October 2015 eruption, OMI also provided a very similar cumulative volume to the field measurements (24.5 × 10 6 m 3 and 23.8 × 10 6 m 3 , respectively). The total lava volume acquired based on OMPS SO 2 emissions was higher (35.3 × 10 6 m 3 ). We note that the approach using fixed c MI S generally gave much lower values (Table 4, Figure 10a,b). Hence, the variation in the pre-eruptive sulfur content estimated above (Figure 9) appears valid, and fixing c MI S at 1050 ppm may be an unrealistic assumption leading to underestimated lava volumes. This is confirmed by looking at the similar estimated total lava volumes for the May 2015 eruption (Table 4). However, we caution that there were UV satellite data gaps during the May 2015 eruption, suggesting a small underestimation for this eruption, as no lava volume was calculated for those days. Remote Sens. 2022, 13, x FOR PEER REVIEW 17 of 26 (a) (b) Figure 10. (a) DRE cumulative lava volume estimated by fixing the sulfur content within melt inclusions at 1050 ppm (35% uncertainty); (b) DRE cumulative lava volume estimated using sulfur contents calculated from the MIROVA lava volume (45% uncertainty). Hence, note the relationship between them. Anomalous sulfur contents (greater than 2000 ppm) were replaced with 1250 or 1600 ppm. As no estimated sulfur content within the melt inclusion was realistic for the April 2007 eruption by using the MIROVA lava volume, the second approach was not applied.

SO2 Flux Measurements from Space and from the Ground-Based NOVAC Network
Comparison between SO2 mass estimation from OMI/OMPS/TROPOMI and SO2 flux from NOVAC measurements illustrates the challenges of comparing satellite-derived masses with ground-based fluxes. The wind speed for the ground-based fluxes, and plume altitude for the satellites, adds significant uncertainty to the values. In addition the vents located on the eastern volcano flank (e.g., February 2019 and April 2020) may be hidden from the NOVAC field of view (Figure 1) when the plume is drifted seawards, hence some SO2 (or the entire plume) may not be observed by the instrument. Plume heights at Piton de la Fournaise range from 1000 m to 5000 m a.s.l. with an average between 3000 m and 3500 m a.s.l. Nevertheless, we find good agreement for the eruptions located within the range of wind directions covered by the NOVAC network (Figures 1  and 3). Adjusting the location of the NOVAC stations and/or expanding the network would be one way to improve plume visibility. Based on our results, moving one or two Figure 10. (a) DRE cumulative lava volume estimated by fixing the sulfur content within melt inclusions at 1050 ppm (35% uncertainty); (b) DRE cumulative lava volume estimated using sulfur contents calculated from the MIROVA lava volume (45% uncertainty). Hence, note the relationship between them. Anomalous sulfur contents (greater than 2000 ppm) were replaced with 1250 or 1600 ppm. As no estimated sulfur content within the melt inclusion was realistic for the April 2007 eruption by using the MIROVA lava volume, the second approach was not applied.

SO 2 Flux Measurements from Space and from the Ground-Based NOVAC Network
Comparison between SO 2 mass estimation from OMI/OMPS/TROPOMI and SO 2 flux from NOVAC measurements illustrates the challenges of comparing satellite-derived masses with ground-based fluxes. The wind speed for the ground-based fluxes, and plume altitude for the satellites, adds significant uncertainty to the values. In addition the vents located on the eastern volcano flank (e.g., February 2019 and April 2020) may be hidden from the NOVAC field of view (Figure 1) when the plume is drifted seawards, hence some SO 2 (or the entire plume) may not be observed by the instrument. Plume heights at Piton de la Fournaise range from 1000 m to 5000 m a.s.l. with an average between 3000 m and 3500 m a.s.l. Nevertheless, we find good agreement for the eruptions located within the range of wind directions covered by the NOVAC network (Figures 1 and 3). Adjusting the location of the NOVAC stations and/or expanding the network would be one way to improve plume visibility. Based on our results, moving one or two of the three stations slightly to the east, or adding a fourth station on the east flank of Piton de la Fournaise would provide better coverage for future eruption vents in this area. On the other hand, the discrepancy between the UV satellite instruments may result from their different spatial resolution and measurement timing leading to either an under-or overestimation relative to NOVAC (Figure 3). Ground-based measurements of plume altitude could be used to reduce the uncertainty in satellite measurements that assume a fixed SO 2 layer height.

April 2007
The April 2007 eruption was characterized by particularly high TADRs and SO 2 fluxes (Figures 4 and A1). Assuming that the lava flow bulk field volume estimated at 185 × 10 6 m 3 by Derrien (2019) [21] is correct, then we see from Table 2 that the volume estimates based on MODIS thermal anomalies generally underestimate (by 32-73 %) the actual volume. As the eruption was particularly intense (peak TADR of 100 m 3 s -1 ), MODIS sensor saturation likely resulted in under-estimated values, meaning that values are "minimum-bounds". This is supported by hot spot temperatures of up to 780 • C [39] and the pixel saturation assessment of Figure 2. Indeed, 1 % coverage by surfaces at 780 • C in a 1 km 2 pixel would have been sufficient to saturate Band 21, equivalent to a 9800 m 2 hot spot. Moreover, some of the lava flows were tube-fed and entered the ocean [39]. Consequently, part of the thermal emission is missing. In addition the ash plume that formed during the peak phase (6-7 April) would have prevented the detection of the hot spots, or dampened the signal. To counter this effect, TADR estimations were increased by 30%. However, this still results in an underestimate. To square with the field-based measurements, the results of Table 2 shows that the adjustment factor needs to be 73%.
On the other hand, using SO 2 emissions we found a larger DRE lava volume (369.8 ± 178.9 × 10 6 m 3 ) than that measured in the field (120.3 × 10 6 m 3 ), suggesting an S excess problem. This excess sulfur has been explained either by the contribution of an unerupted amount of deep S-rich magmas [13] or by release from the hydrothermal system. However, according to our results using seismically-derived TADR data [34], we estimated an average S within melt inclusion of between 1050 and 1250 ppm (Figure 7a), corresponding to petrological analysis [36]. Hence, we can conclude that no S excess was associated with this eruption and that suggestions of an excess arose from preliminary underestimates of erupted lava volume due to saturation issues. This paroxysmal eruption released a large amount of SO 2 ( Figure A1) due to the fast ascent of very deep magma from the crustal underplating depth [13,34]. This deep magma input would have pressurized the shallow reservoir and led to this type of intense activity, as suggested by Di Muro et al. (2014) [36]. Following Walker (1988) [40], that the summit caldera collapses were preceded by an increase in TADR suggests an enhancement in the drainage of the shallow system, leading to a passive form of caldera collapse. This behavior is observed at other basaltic shields that experience high effusion rate and/or voluminous eruptions (e.g., Kilauea (Hawaii [41,42])).

May 2015
This eruption started on May 17 and was characterized by high initial TADR (30 m 3 s -1 ) followed by an exponential decay to 14 m 3 s -1 (Figure 4). The same trend is evident in the SO 2 emissions (Figures A1 and 4a) and is typical of tapping of a pressurized source [43]. We also point out that this eruption has a similar TADR trend as the April 2007 eruption ( Figure 5). Both eruptions show an exponential decay of TADR despite different eruptive processes. In fact, the May 2015 eruption was fed by a closed system, and that observed in the much larger 2007 eruption was fed by an open system. This poses the possibility of using TADR trends to constrain deep magma transfer ( Figure 5). According to field measurements, the bulk volume erupted in May 2015 was 5.7 × 10 6 m 3 . The total lava volume derived from MODIS data analysis is in good agreement except for the MODVOLC data, where we obtained a lava volume of 9.0 ± 2.7 × 10 6 m 3 . Despite MODVOLC data gaps, this system overestimates the total lava erupted, which we assume is mostly due to the bow-tie effect (Figure 11). At high scan angles, Earth's curvature cannot be neglected by the MODIS instrument, leading to scan-to-scan overlap so that hot spot radiances are counted twice [44]. lava volume derived from MODIS data analysis is in good agreement ex MODVOLC data, where we obtained a lava volume of 9.0 ± 2.7 × 10 MODVOLC data gaps, this system overestimates the total lava erupted, whic is mostly due to the bow-tie effect (Figure 11). At high scan angles, Earth's cu not be neglected by the MODIS instrument, leading to scan-to-scan overlap spot radiances are counted twice [44]. Using OMI and OMPS SO2 data, the total DRE volume estimated using eruptive S content of 1050 ppm underestimates the lava emitted giv × 10 m 3 (Table 4). However, using a lower of 430 ± 157 ppm yields vo closer to that measured in the field (3.5 × 10 m 3 ) (Figure 10b). Note that underestimation may also be an artifact of the low temporal resolution of the and/or an inaccurate assumption of SO2 plume altitude, leading to an unre time-series of SO2 emissions. According to the reverse petrological analysis, t tive sulfur content of 430 ± 157 ppm may be slightly too low (Figure 8). A eruptive sulfur contents of 500 ppm are known, such cases are rare [33]. How chemical studies of the evolved and highly crystallized 2014 magma have sho of the pre-eruptive S content can be locked in solid phases (sulphides) and no and that can further influence the amount of sulfur immediately available to t upon magma ascent and decompression [45].
We also note that the May 2015 eruption emitted a more evolved magm [46]. In mid-April 2015, deep seismicity and an increase in CO2 discharge we suggesting the ascent of new mafic magma from depth into the shallow rese 1.5 km below sea level [47]. Petrological studies show that magma cooling an were effective since 2014 or earlier and that shallow seated May 2015 magmas the most evolved and lacked any evidence of interaction with hotter and d [43]. No apparent mafic input is recorded in the sulfur content variation in F gesting that magma supplying the May 2015 eruption came from the uppe shallow plumbing system. This was extruded due to the arrival of new magm plumbing system but without interacting with it. This could explain the low constant sulfur content implied by the reverse petrological approach. This is a mon process at basaltic shields [48]. Using OMI and OMPS SO 2 data, the total DRE volume estimated using a fixed preeruptive S content c MI S of 1050 ppm underestimates the lava emitted giving~0.6-1.1 × 10 6 m 3 (Table 4). However, using a lower c MI S of 430 ± 157 ppm yields volumes much closer to that measured in the field (3.5 × 10 6 m 3 ) (Figure 10b). Note that any volume underestimation may also be an artifact of the low temporal resolution of the UV sensors and/or an inaccurate assumption of SO 2 plume altitude, leading to an unrepresentative time-series of SO 2 emissions. According to the reverse petrological analysis, the pre-eruptive sulfur content of 430 ± 157 ppm may be slightly too low (Figure 8). Although pre-eruptive sulfur contents of 500 ppm are known, such cases are rare [33]. However, petrochemical studies of the evolved and highly crystallized 2014 magma have shown that part of the pre-eruptive S content can be locked in solid phases (sulphides) and not in the melt and that can further influence the amount of sulfur immediately available to the gas phase upon magma ascent and decompression [45].

August-October 2015
We also note that the May 2015 eruption emitted a more evolved magma than usual [46]. In mid-April 2015, deep seismicity and an increase in CO 2 discharge were recorded, suggesting the ascent of new mafic magma from depth into the shallow reservoir at 0.5-1.5 km below sea level [47]. Petrological studies show that magma cooling and degassing were effective since 2014 or earlier and that shallow seated May 2015 magmas were among the most evolved and lacked any evidence of interaction with hotter and deeper inputs [43]. No apparent mafic input is recorded in the sulfur content variation in Figure 9, suggesting that magma supplying the May 2015 eruption came from the upper part of the shallow plumbing system. This was extruded due to the arrival of new magma inside the plumbing system but without interacting with it. This could explain the low and nearly constant sulfur content implied by the reverse petrological approach. This is another common process at basaltic shields [48].

August-October 2015
The August-October 2015 eruption can be divided into three phases on the basis of TADR (Figure 4a). The first phase was marked by a high TADR (peak of 50 m 3 s -1 ) and was followed by a rapid decline to a nearly constant TADR at 5 m 3 s -1 in phase 2. Three short-lived pulses, each about two days long, characterized the last phase. The bulk lava volume erupted during this last 6-day long phase represented almost 50 % (∼11 × 10 6 m 3 ) of the total emitted from August 24 to October 17 (∼32 × 10 6 m 3 ). By integrating the TADR through time, the MIROVA and Manual methods yield equivalent total volumes. However, the MODVOLC system again overestimates the volume by~40 × 10 6 m 3 ( Table 2). As with the May 2015 eruption, we obtain an underestimation by 13 to 20 × 10 6 m 3 of the total DRE lava volume if we fix the pre-eruptive sulfur content at 1050 ppm (Table 4). Depending on the UV sensor, we obtained 3.9 ± 1.8 × 10 6 m 3 (OMI) and 10.5 ± 4.8 × 10 6 m 3 (OMPS) when we expect a total lava volume of 23.8 × 10 6 m 3 according to field measurements. Using the pre-eruptive sulfur content from the reverse petrological approach (± 466 ppm), DRE volume estimations are better (OMI: 24.5 ± 17.2 × 10 6 m 3 ). The OMPS dataset slightly overestimates (35.3 ± 24.7 × 10 6 m 3 ) the lava emitted, likely due to its lower spatial resolution. Figure 9 reveals a nearly constant c MI S of 200-700 ppm before an increase to 2000 ppm during the last phase. This may suggest an evolution from a degassed melt during the two first phases to a mafic gas-rich melt during the last stage of the event.
Because of the involvement of highly degassed evolved melts, the initial phase of strong lava flux does not correspond to the most intense phase of gas emissions (Figures 4 and A1). This is consistent with Sundermeyer et al. (2019) [46], who modelled diffusion times within olivine crystals. Sundermeyer et al. (2019) [46] noted a progressive increase of the MgO magma content during the third phase, which was accompanied by increased CO 2 in summit fumarole emissions [12]. In addition, the bulk composition of erupted products indicate that recharging magmas dominated evolved melts during the final stage, probably because the shallow reservoir had been almost emptied by the three previous eruptions in 2015. That could be explained by the "Stromboli effect" [12]. According to , after emptying the shallow part of the plumbing system, the deeper part arrives at the surface producing an increase in TADR and SO 2 emissions [12].

February 2019
The trends in TADR and SO 2 emissions (Figure 4) show this eruption to be characterized by a terminal burst. This was related to the opening of a new E-W trending fissure on 7 March, 17 days after the onset of the eruption and 3 days before the end. TADR was estimated at~10 m 3 s -1 before this event increasing to 19-37 m 3 s -1 on 9-10 March (Figure 4a). A peak in SO 2 emissions is observable with OMI (4.9 kt) and OMPS (19.2 kt) on 9 March, coinciding with the TADR increase ( Figure A1). According to the bulk field data, the total lava volume emitted during this event was around 14.5 × 10 6 m 3 . The MIROVA dataset provides a similar value (Table 2). Despite the absence of hot spots after 3 March, missing the intense surface activity at the end of the eruption, MODVOLC also gives a reasonable cumulative lava volume estimate (11.2 × 10 6 m 3 ), probably due to the bow-tie effect [44]. On the other hand, the estimate from manual MODIS processing is lower (9.5 × 10 6 m 3 ), meaning that the conversion coefficients of Equation (5) need to be slightly adjusted for this style of activity. Figure 10 shows two different trends in cumulative volume acquired from the reverse petrological approach. Using a fixed c MI S , we observe lower DRE lava volume estimations. However, by varying the pre-eruptive sulfur content, we obtain cumulative lava volumes closer to the field-based measurements. This suggests that a relatively constant c MI S of 528 ± 284 ppm during the eruption is reasonable. As for the 2015 eruptions, this value would be consistent with a degassed differentiated melt coming from the shallow reservoir. However, we observe an increase in apparent c MI S at the end of the eruption, corresponding to the increase in TADR and SO 2 emissions. This may suggest that the fissure eruption was fed by a volatile-rich, less evolved melt: formation of a new dyke or emptying of the shallow reservoir allowing deep undegassed magma to erupt as in the third phase of the August-October 2015 eruption, as suggested by   [12].

April 2020
The 2-6 April 2020 eruption occurred close to the site of the February 2019 eruption. Weak ground deformation and a low number of seismic events suggest that the pathway for magma ascent was already open following the 2019 eruption [35]. On 4 April, seismicity increased until 6 April, when the eruption stopped abruptly. During this period, large quantities of Pele's hair were emitted [35], indicating a unusual highly explosive event for Piton de la Fournaise feeding high fountains and associated with strong volcanic degassing. This is supported by the total mass of SO 2 detected during 4-6 April by TROPOMI, OMPS, and OMI (17.32 kt, 36.67 kt, and 4.12 kt, respectively). An OMI data gap on 6 April during the highest SO 2 emissions ( Figure A1) explains the low OMI value. According to Peltier et al. (2020) [35], the bulk lava volume ranges from 6-10 × 10 6 m 3 . MODIS image processing yields a bulk volume of 2-4 × 10 6 m 3 with a~35% error ( Table 2). Cloud cover significantly compromised the availability of usable (cloud-free) images, indicating that caution should be used when using satellite-derived data, which should always be cloud-screened to check for cloud-cover induced trends and drop-outs [49].
On the other hand, the DRE lava volumes estimated using the reverse petrological approach are higher ( Table 4). The cause of the high bias in the OMPS-derived volume (nearly 20 × 10 6 m 3 ) is unclear and requires further investigation. In addition, the results indicate that fixing c MI S at 1050 ppm is a reasonable assumption for this eruption. According to the petrological data ( [36]; Figure 8), a sulfur content of 1050 ppm represents, on average, the sea level magma reservoir of Piton de la Fournaise, suggesting that the April 2020 eruption was supplied by this main reservoir. However, Figure 9 also shows an increase in pre-eruptive sulfur content, probably resulting from deeper mafic inputs. The continuous increase in soil CO 2 fluxes and syn-eruptive seismicity supports the idea of a deep magma influx into the sea level reservoir [35] after initial emptying.

Conclusions
Multiple datasets were acquired for five eruptions at Piton de la Fournaise, La Réunion. Analysis of the temporal evolution of TADR, combined with the associated SO 2 emissions and inferred pre-eruptive magmatic sulfur content, reveals that Piton de la Fournaise eruptions may follow several distinct trends, including classic exponential decay (April 2007, May 2015) and terminal burst (August-October 2015, February 2019 and April 2020). In addition, our study confirms no evidence of volatile excess associated with the eruptions considered, which we conclude were mostly fed by shallow degassed magma; instead, we identify a volatile deficit in most cases.
Manual processing of MODIS data validates the efficiency of hot spot detection and TADR-derivation by the MIROVA system during the effusive eruptions at Piton de la Fournaise, meaning that the conversion coefficient in Equation (5) is valid. In contrast, we find that the MODVOLC system often overestimates TADR and also the total DRE erupted volume. This is unexpected given that the MODVOLC system often fails to detect hot spots, which should lead to underestimation of TADRs and lava volumes. This is particularly apparent for the April 2020 eruption. For this case, despite having only three measurements (compared to 11 for MIROVA and manual processing), MODVOLC provides the highest DRE lava volume estimates, which is likely due to double counting resulting from the MODIS bow-tie effect. Nevertheless, lava erupted during bad weather is significantly underestimated by all the approaches compared to the field measurements. Underestimation of lava volume also occurred during the paroxysmal eruption in April 2007 due to widespread saturation of the MODIS TIR channels ( Table 2).
Using the reverse petrological method and the erupted lava masses obtained from MODIS data, we derived time-series of pre-eruptive sulfur contents, which were compared with sulfur concentrations measured in melt inclusions. Some results show lower values than expected from typical undegassed magmas. This could be explained by the involvement of shallow seated, evolved, and partly degassed magma in the eruptions. On the other hand, increased TADR and SO 2 emissions at the end of the February 2019 and August-October 2015 eruptions (terminal bursts) suggest an increase in pre-eruptive sulfur content. This may indicate replenishment of the shallow magma reservoir by deeper, volatile-rich mafic inputs. In addition, total DRE lava volumes estimated using the SO 2 emissions detected by OMI, OMPS, and TROPOMI are in good agreement with ground-truth. Hence, using the SO 2 emissions and pre-eruptive sulfur contents could be a viable alternative to estimate bulk lava volumes during bad weather, when the satellite-derived TADRs may be compromised. However, we also recognize that during very bad weather conditions when both thermal anomalies and SO 2 plumes are partially or wholly undetectable, neither satellite-based technique (IR or UV) may provide accurate results.
The analysis of multiple datasets for a sequence of effusive eruptions from a single magmatic system at Piton de la Fournaise allowed the identification of two groups of eruptive behavior. The first group is characterized by a strong initial peak in TADR and SO 2 emissions before an exponential decrease for the classic "waxing-waning" trend for effusive eruptions defined by Wadge (1981) [43]. Note that the intensity of the SO 2 peak will depend on the degree of pre-eruptive degassing of the magma. This trend is represented by two eruptions driven by different eruptive processes (April 2007: piston-collapse [34] or deep input [33] or a combination of both; May 2015: mostly volatile exsolution). Hence, it would be interesting to consider in detail the different sources of overpressure leading to this eruptive behavior. The second group (the August-October 2015, February 2019 and April 2020 eruptions) is characterized by an intense and explosive terminal phase. For these cases, we show that the late increase in TADR (terminal burst) is most likely related to the arrival of new, deep undegassed magma causing an increase in gas emissions (e.g., enrichment in CO 2 within the summit fumaroles, soil, and SO 2 emissions). However, the driving process for this type of behavior remains unclear and is yet to be clearly identified.
Based on this work, we suggest three improvements to aid with the remote sensingbased monitoring of Piton de la Fournaise and some potential areas of future research: 1.
The addition of another NOVAC DOAS station on the east flank of Piton de la Fournaise will allow efficient surveys of future vents in this area. In addition, a more extended comparison between the satellite SO 2 masses and the ground-based SO 2 fluxes could allow cross-validation of both methods, depending on measurement conditions. A combination of both would lead to better results.

2.
This study demonstrates that the MODVOLC did not detect a few hot spots in addition to overestimating the TADR and hence the total lava volume. Consequently, MOD-VOLC data should be used with caution when used for TADR conversion. However, MIROVA appears well calibrated to produce reliable TADR for Piton de la Fournaise.

3.
In addition, depending on the type of activity (e.g., purely effusive, fountains, etc.), the residual sulfur content in the matrix can vary from 50 ppm up to 230 ppm. Hence, our use of a fixed value of 160 ppm for the development of time-series of pre-eruptive sulfur contents is a significant assumption. Future work may select the optimal value according to the style of activity.
Finally, we also suggest further analysis of ground deformation and seismicity datasets or use of Bayesian inversion methods to model the source of deformation during magma migration at Piton de la Fournaise [50]. These data could improve our understanding of the different eruptive trends observed in this study. Funding: S.A. is thankful for support from the Swedish National Space Agency, grant Dr. 149/18. S.C. acknowledges NASA support from the Aura Science Team (grants 80NSSC17K0240 and 80NSSC20K0983) and the Science of Terra, Aqua and Suomi/NPP program (grant 80NSSC18K0688). This is ANR-LAVA publication number 21, and ClerVolc publication number 519.

Data Availability Statement:
The data that support the findings of this study are available from the authors upon reasonable request.