Francisco J Vasconez

Silvana Hidalgo

Jean Battaglia

Stephen Hernandez

Patricio Ramón

Benjamin Bernard

Diego Coppola

Sébastien Valade

Santiago Arellano

Céline Liorzou

Marco Almeida

Marcelo Ortíz

Jorge Córdova

Anais Vásconez Müller

Linking ground-based data and satellite monitoring to understand the last two decades of eruptive activity at Sangay volcano, Ecuador

Keywords: Sangay volcano, long-lasting activity, episodic eruptions, size of the eruptions, eruption frequency

INTRODUCTION

Sangay (78.34°W; 2.00°S; 5326 m asl) is the southernmost active volcano of the Northern Volcanic Zone (NVZ) in the Andes (Fig. 1a). It lies in a remote area of the eastern flank of the Cordillera Real, within a rainforest with no car access and where rainy and cold weather conditions are prevalent during most of the year. Therefore, accessing the volcano is logistically difficult, giving way to only few direct observations. Sangay is a conical stratovolcano, which stands out by having a dome and three main craters at the summit area and a lateral vent, aligned west-east (Fig. 1b). According to [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF] these vents evolved through time as follows: in 1628 A.D. there was only one Central Crater (CC) through which lava flows and eruptive plumes were emitted (Fig. 1b). This crater remains active until today. Later, in 1956, a new crater located to the Northeast (NEC) of the Central Crater was reported, and this remained active until 1976. During 1974 and1975 two new vents also appeared at the western zone of the Central Crater: the West Crater (WC) and West dome (W-dome). The latter began as a fissure and became a lava dome in 1983. Activity of this dome claimed the life of two English explorers back in 1976 (Global Volcanism Program 1976). Finally, the Ñuñurcu vent (ÑV) is a lateral vent located to the southeast of the Central Crater from where in 2007 a lava dome/coulée was extruded. This vent was active until June 2020 and extruded lava flows down the southeastern flank of the edifice. Moreover, two large landslide scars account for the main morphology of Sangay volcano [START_REF] Valverde | Enormous and far-reaching debris avalanche deposits from Sangay volcano (Ecuador): Multidisciplinary study and modeling the 30 ka sector collapse[END_REF]. These scars and their associated debris avalanches occurred between 250-100 ka BP and 29 ka BP, respectively [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF][START_REF] Valverde | Enormous and far-reaching debris avalanche deposits from Sangay volcano (Ecuador): Multidisciplinary study and modeling the 30 ka sector collapse[END_REF], and both affected the eastern and southeastern flanks leaving debris avalanche deposits where the city of Macas is now located with more than 24,400 inhabitants (Villareal 2017; Valverde et al. 2021, Fig. 1c).

Fig. 1 a) Main geodynamic features of Ecuador, which include Galápagos Spreading Center (GSC) and GFZ (Grijalva Fracture Zone); the red star shows the location of Sangay volcano, while Chimborazo and Tungurahua volcanoes and the main cities are displayed as the orange triangles and the black dots, respectively. b) Aerial view of Sangay with its main features: West-dome (W-dome), West Crater (WC), Central Crater (CC), Ñuñurcu vent (ÑV) and Northeast Crater (NEC). The ages of the active vents are mainly based on [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF]. Chimborazo volcano is visible in the background. c) Location of the 10 broadband seismic stations used to detect Sangay seismicity, identifies with 4-letter codes. SAGA station and the regional seismic network are pointed out as colored dots, while main rivers and cities are highlighted in blue and black, respectively. The location of Tungurahua, Chimborazo and Sangay volcanoes are also shown.

Since 1628 A.D. eruptive activity at Sangay has been mostly of continuous nature (Hall 1977). The Global Volcanism Program (1996) reports eruptions from 1728 A.D. to 1916 A.D. and again since 1934 A.D., seemingly more or less continuous. Based on historical records, in 1628, Riobamba city, located 50 km to the northwest of the volcano (Fig. 1c), was completely blanketed and darkened by ash for one and a half days, disrupting the daily life of the society of that epoch [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF]. Additionally, over the last two decades and despite the almost continuous activity, on a few occasions has ash fallout affected large cities such as Guayaquil in 2018 and 2020, 170 km away from Sangay [START_REF] Igepn | Special report N o 2, Sangay: Continua el proceso eruptivo[END_REF][START_REF] Igepn | Special report N o 3, Sangay: Actualización de la actividad eruptiva[END_REF]. The low frequency of ash fallout reaching major cities is mainly due to the small height of the eruptive plumes, the low amount of ash and the remote location of the volcano. However, small communities sparsely located ~30-40 km to the west of the volcano, which is the prevalent direction of the winds [START_REF] Bernard | Dispersion of volcanic ash clouds in Ecuador: a 20 years perspective[END_REF], have been affected by ash fallouts leading to the loss of crops and causing cattle diseases. Potential aviation hazards also exist given that the eruptive plumes and co-pyroclastic ash plumes could traverse air traffic routes near Sangay. Other phenomena like lava flows, pyroclastic density currents, and lahars are also common, but they are usually restricted to the volcanic ring plain (< 10 km away from the summit). All these phenomena were synthesized on the first hazard map published by [START_REF] Ordoñez | Multi-year regional infrasound detection of Tungurahua, El Reventador, and Sangay volcanoes in Ecuador from 2006 to[END_REF]. Finally, debris avalanches are possible to occur as Sangay grows higher and more unstable over time [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF].

Throughout the years, only few scientific studies have been carried out on Sangay. [START_REF] Johnson | Plugs and chugs-seismic and acoustic observations of degassing explosions at Karymsky, Russia and Sangay, Ecuador[END_REF], [START_REF] Konstantinou | Nonlinear Time Series Analysis of Volcanic Tremor Events Recorded at Sangay Volcano, Ecuador[END_REF], and [START_REF] Lees | Non-linear explosion tremor at Sangay, Volcano, Ecuador[END_REF] Furthermore, [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF] performed an impressive field mission to collect samples from this remote volcano constructed over the last 500 ky. Interestingly, Sangay displays a particular geochemical signature when compared with the rest of the volcanoes in the Ecuadorian arc. They found Niobium-rich rocks and the most primitive basalt. In addition, the olivine-hosted melt inclusions analyzed from that basalt show high volatile (Sulfur-Chlorine-Fluorine) contents pleading for a volatile-enriched source for the magma [START_REF] Narvaez | Constraining magma sources using primitive olivine-hosted melt inclusions from Puñalica and Sangay volcanoes (Ecuador)[END_REF]. In a broader context, the magmatism in the southern termination of the Ecuadorian volcanic arc is the result of the partial melting of the mantle metasomatized by aqueous fluids released by an old and cold slab beneath the continental South American plate (Fig. 1a, [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF][START_REF] Narvaez | Constraining magma sources using primitive olivine-hosted melt inclusions from Puñalica and Sangay volcanoes (Ecuador)[END_REF].

In this study, we used data from 10 permanent seismic stations, one of them coupled to an acoustic sensor and a ScanDOAS (Differential Optical Absorption Spectroscopy, Fig. 1c), instrument from the Network for Observation of Volcanic and Atmospheric Change (NOVAC). All these instruments were installed and are operated by the Instituto Geofísico (IG-EPN). This information is complemented by satellite observations to monitor the most recent eruptive activity of this remote volcano. The purpose of the present contribution is to address the lack of knowledge about Sangay's volcanic activity by linking the available ground-based geophysical data and the satellite imagery of the last two decades. These data allow us to characterize internal and superficial activity as well as to quantify, within the intrinsic methods' uncertainties, the size of the eruptions in terms of magma output rates and cumulative volumes based on thermal satellite information. We also hypothesize some mechanisms to explain Sangay's long-term eruptive activity based on the extruded lava volume, the variations of the whole-rock geochemistry of the tephra products and the recorded geophysical parameters.

METHODS AND DATA PROCESSING

Ground-based monitoring network

Given the low risk associated with Sangay's activity and the extremely difficult access, the permanent monitoring network is modest. A single broadband seismometer (SAGA-BB) coupled with an acoustic sensor and a DOAS scanning instrument were installed 6 km from the summit in September 2013 (Fig. 1c).

Nevertheless, extreme weather conditions and frequent ash fallout make it very difficult to keep these instruments working during long periods of time. As a result, the SAGA-BB has operated intermittently from September 2013 to the present, with no data during periods as long as a few months. Furthermore, the seismic digitizer had to be replaced in November 2018 due to complete damage. Consequently, most of the seismic activity at Sangay had to be complemented by nine regional broadband seismic stations located from 56 to 102 km away from the volcano (Fig. 1c). On the other hand, the DOAS scanning instrument provided good quality information only between December 2013 and June 2015 and again since December 2019 when the spectrometer was replaced.

In order to quantify the seismicity at Sangay, we calculated several parameters at SAGA-BB station and at the nine regional seismometers (Fig. 1c). For SAGA-BB, we first applied an STA/LTA (short-term average amplitude/long-term average amplitude) algorithm to count the number of transient events following the approach presented by [START_REF] Battaglia | Families of similar events and modes of oscillation of the conduit at Yasur volcano (Vanuatu)[END_REF]. We used a STA of 1 second and a LTA of 60 s, applied to data filtered between 0.5 and 20 Hz. Peak-to-peak amplitudes were calculated on-the-fly for each detection in the same frequency band. Secondly, to identify families of similar events, we extracted 20.48s signal windows for the vertical component of SAGA-BB, for all detections with an STA/LTA higher than 3.0. In total, we compared nearly 64,000 waveforms using cross-correlation, after filtering the waveforms between 0.8 and 30 Hz. To build families we used open clustering allowing chain similarity [START_REF] Battaglia | Families of similar events and modes of oscillation of the conduit at Yasur volcano (Vanuatu)[END_REF] and a correlation threshold of 0.80.

To complement and fill in observational gaps during SAGA outages, we also designed and implemented a subspace detector using the vertical component of the nine regional broadband seismometers (Fig. 1c).

Subspace detectors are a class of correlation detectors but with the added benefit of being less sensitive to subtle changes in source-time history, mechanism, or location [START_REF] Barrett | An empirical approach to subspace detection[END_REF][START_REF] De | Empirical subspace detection applied to triggered seismicity by the July 25, 2011, Mw 5.0 earthquake in the Sea of Marmara, Turkey[END_REF]. We employ the strategy of [START_REF] Harris | Subspace detectors: theory[END_REF] to generate a design matrix of suitable events emanating from the volcano. For each regional sensor we filter, align, and build our basis functions (templates) from these matrices via singular value decomposition. Because of the large distances between this regional network and the volcano, and the generally low magnitude nature of the explosions, the signal is only coherent across a very narrow frequency range: we find that the signal to noise ratio is maximum in the 0.6 -1.2 Hz range. Our detector calculates a sufficient statistic (akin to a correlation coefficient) using an empirically-determined total of five basis functions. A detection is declared when the sufficient statistic surpasses a threshold of 0.1 on at least three channels. The selection of 5 basis functions and a detection threshold of 0.1 represents the best tradeoff between code execution time and rate of false positives. This procedure was run from August 2007 to May 2020 generating a total of ~238,000 detections.

Unfortunately, because of the extremely narrow nature of the filter centered over 1 Hz, raw detections from the subspace detector are contaminated with numerous false detections from not only random noise but also local, regional, and even teleseismic earthquakes. To address this problem, we applied a binary classification algorithm written in MATLAB and using gentle adaptive boosting. This classifier was trained and validated on a suitably large dataset of true positives and known spurious detections (false positives).

The training and validation sets were mutually exclusive to prevent data leakage, with respective labels being confirmed using the nearfield sensor SAGA-BB as ground truth. A total of 119 features for each event were extracted from the waveforms themselves and, consisted of various standard statistical measures such as mean, rms, kurtosis, amongst others. In the end, the most important features were the sufficient statistic itself, and the cross-correlation lag values (between sensors) of candidate events. After validation, we find that the classifier has a true positive identification rate of >97%, and false positive identification rate of >96%, which we find suitable for the purposes of this study since we seek only broad event detection capability for those periods of time when SAGA-BB is not transmitting. Finally, we apply the classifier to the entire dataset of ~238,000 detections, and we were able to prune this number down to ~98,000 total events detected between 2007 -2020. This procedure constitutes a significant improvement in the homogeneity of the Sangay seismic catalog, albeit at the cost of a larger magnitude of completeness; in general we find that the subspace detector coupled with the binary classifier can accurately identify the majority of Sangay events with amplitudes greater than ~8 m/s as recorded on the SAGA instrument, with anything lower than 8 m/s being too low in amplitude to properly emerge from the noise (Supplementary Material 1).

Finally, SO2 daily fluxes were calculated using the NOVAC software [START_REF] Johansson | Application of passive DOAS for studies of megacity air pollution and volcanic gas emissions[END_REF] introducing the wind conditions provided by NOAA wind models (National Oceanic and Atmospheric Administration, https://www.ready.noaa.gov/READYamet.php) calculated at the summit altitude of the volcano, assuming that the plume drifts horizontally at the vent altitude. Whenever Volcanic Ash Advisory Center (VAAC) alerts were available, we used the reported plume height, plume direction and velocity instead of the NOAA data. In addition, we calculated daily observed SO2 emission following the approach of [START_REF] Hidalgo | SO2 degassing at Tungurahua volcano (Ecuador) between 2007 and 2013: Transition from continuous to episodic activity[END_REF]. This approach is justified given that the dominant wind direction is towards the station.

Satellite observations and overflights

Satellite-based instruments are very useful tools for monitoring active volcanoes, in particular for those located in remote zones, as is the case of Sangay. We used the freely available data provided by the following systems: (1) Middle InfraRed Observation of Volcanic Activity (MIROVA, https://www.mirovaweb.it/?action=volcanoes) and Fire Information for Resource Management System (FIRMS, https://firms.modaps.eosdis.nasa.gov) that provided thermal information, (2) Washington VAAC alerts, which supplied information about the dispersal and height of ash plumes (https://www.ssd.noaa.gov/VAAC/archive.html), (3) SO2 mass derived from TROPOMI (TROPOspheric Monitoring Instrument) processed by Monitoring Unrest from Space (MOUNTS, http://www.mountsproject.com/home), and from OMI sensor processed by National Aeronautics and Space Administration (NASA) via Global Sulfur Dioxide Monitoring (https://so2.gsfc.nasa.gov/), and (4) overflights to capture the volcanic products, in particular, lava flows and pyroclastic density currents, whenever the weather conditions allow to obtain clear images.

We used an automated global hot spot detection system based on real-time Moderate resolution Imaging Spectroradiometer infrared data (MODIS), which then is processed by MIROVA [START_REF] Coppola | Lava discharge rate and effusive pattern at Piton de la Fournaise from MODIS data[END_REF][START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF]. MIROVA completes automatic detection and location of thermal anomalies and provides a quantification of the Volcanic Radiative Power (VRP) in watts (W) [START_REF] Coppola | Lava discharge rate and effusive pattern at Piton de la Fournaise from MODIS data[END_REF][START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF][START_REF] Coppola | Enhanced volcanic hot-spot detection using MODIS IR data: results from the MIROVA system[END_REF][START_REF] Coppola | Modelling satellite-derived magma discharge to explain caldera collapse[END_REF].

From VRP we were able to derive Volcanic Radiative Energy (VRE), time-averaged discharge rate (TADR) and the cumulative lava volume according to the methodology proposed by several authors [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF][START_REF] Coppola | Modelling satellite-derived magma discharge to explain caldera collapse[END_REF][START_REF] Naismith | Eruption frequency patterns through time for the current (1999-2018) activity cycle at Volcán de Fuego derived from remote sensing data: Evidence for an accelerating cycle of explosive paroxysms and potential implications of eruptive activity[END_REF][START_REF] Walter | Complex hazard cascade culminating in the Anak Krakatau sector collapse[END_REF]. Volcanic Radiative Energy is the energy in Joules obtained by the trapezoidal integration of the Volcanic Radiative Power timeseries. Weekly time-averaged discharge rate is obtained from the average VRP during a period of seven days and the chemical composition of the extruded lavas. For this, we used the silica content of samples collected and analyzed in the 1990s [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF], and samples from this study corresponding to eruptive phases of 2015, 2018, 2019 and 2020. Finally, the cumulative volume was estimated by using the weekly time-averaged discharge rate and a week time-integration, which we considered as the best approach at volcanoes that have strong cloud cover, as is the case of Sangay.

Washington VAAC reports alerts for Sangay since January 2004. Since then, 900 confirmed alerts have been reported until 31 May 2020. Moreover, TROPOMI sensor [START_REF] Theys | Global monitoring of volcanic SO2 degassing with unprecedented resolution from TROPOMI onboard Sentinel-5[END_REF], Hedelt et al. 2019) on board Sentinel-5P provides SO2 vertical column densities maps since 2018. The volcano monitoring system MOUNTS [START_REF] Valade | Towards Global Volcano Monitoring Using Multisensor Sentinel Missions and Artificial Intelligence: The MOUNTS Monitoring System[END_REF]) converts these into daily SO2 masses, with a spatial filter allowing its attribution to Sangay. Finally, sporadic overflights (one or two per year) using a visual and a FLIR thermal camera complemented with satellite images (Suomi-NPP, Aqua and Terra) provided by the Fire Information for Resource Management System, give information on the emplacement of lavas and pyroclastic density current deposits. In order to filter out wrong satellite location of the thermal anomalies because of the angle between the satellite flight and the topography, we applied a <0.5 threshold on the satellite-track data [START_REF] Wang | Assessment of the NOAA S-NPP VIIRS Geolocation Reprocessing Improvements[END_REF].

Volcanic products geochemistry

We analyzed bulk ash samples from the 2015, 2018, 2019 and 2020 eruptions collected at SAGA-BB site (Fig. 1c). Major and trace elements were determined at the PSO/IUEM (Pôle Spectrométrie Océan, Institut Universitaire Européen de la Mer, Brest, France), following the analytical procedure of [START_REF] Cotten | Origin of anomalous rare-earth element and yttrium enrichments in subaerially exposed basalts: evidence from French Polynesia[END_REF].

Typically, 250 mg of rock power were dissolved in closed screw-top teflon vessels (Savillex) at about 90°C for one day using 3 ml of concentrated HF, and 1 ml of concentrated HNO3. Next, 96 ml of H3BO3 aqueous solution (20 g/L H3BO3) were added to neutralize the excess HF. All reagents used are analytical grade.

Elements were measured by inductively coupled plasma-atomic emission spectrometry (ICP-AES) using a Horiba Jobin Yvon® Ultima 2 spectrometer. The boron included in the solution was used as an internal standard. Calibrations were made using international standard, ACE, WSE, JB2. For major elements, relative standard deviation is ≤1% for SiO2 and ≤2% for the other major elements. Finally, we used for comparison the geochemical data of most recent products (1990s) from Sangay's summit area reported by [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF] (Supplementary Material 2).

RESULTS

Seismic activity

Local seismic observations

Despite the intermittence of SAGA-BB (Fig. 2a), we were able to identify quiescent and eruptive phases based on the number and amplitude of the seismic events recorded from September 2013 to May 2020 (Fig. 2b &2c). Quiescent phases are characterized by a low number of detections, mostly related to background noise or regional events (Fig. 2b). In contrast, eruptive phases display clear peaks above background level in both number of events and amplitude, reaching more than 100 seismic events per day and amplitudes above 60 m/s (Fig. 2b &2c). In addition, we observed that there is no evident escalating unrest before eruption onsets (Fig. 2b). We identified four major seismic families on SAGA-BB which occurred successively and only during eruptive periods (Fig. 2b &2c). They are active only one at a time. The transition between families A and B (August 2019) is likely caused by a change in the instrument response as the station was suddenly tilted and flooded at the time of the transition. On the other hand, the two later transitions (March and April 2020) may be related to changes in the source or morphology of the volcano (Fig. 2d). Regional stations are too distant to record the finer details in the waveforms that SAGA-BB provides, though some broad degradation in waveforms similarity is discernable and match to a degree the observations seen at SAGA-BB, in particular the March 2020 transition (Fig. 3a).

Overall, many of the events belonging to these families display acoustic waveforms on the infrasound channel, indicating that they are explosion quakes. Figure 2d displays monthly stacked waveforms showing a great similarity over years but also a significant change in the waveform and frequency content since March 2020. These events can be interpreted as the resonance of the conduit and surrounding volcanic edifice in response to the explosive decompression of the magma [START_REF] Battaglia | Families of similar events and modes of oscillation of the conduit at Yasur volcano (Vanuatu)[END_REF]. Details of these modes of oscillation are filtered by the distance, inducing the presence of a single family during most of the study period. Some of the explosion quakes are also followed by long oscillations, sometimes monochromatic, possibly corresponding to chugging following the initial explosions whose waveforms are identical to shorter events, similar to what is described by [START_REF] Johnson | Plugs and chugs-seismic and acoustic observations of degassing explosions at Karymsky, Russia and Sangay, Ecuador[END_REF]. The families group between 0 and 90% of the daily larger events (STA/LTA > 6) during eruptive periods with a global average at 4024%. This indicates the presence of a significant and dominating stable explosive process during these phases. A general tendency to decrease toward the end of the study period is observed. Further investigation is needed to better understand the changes that occurred in March 2020.

Regional seismic network

Near-field seismic monitoring of Sangay volcano has proven difficult to keep running continuously due to its remote location and frequent cloud and ash coverage, although efforts are currently underway to expand the network and improve the reliability of existing instrumentation. In the meantime, the subspace detector described above has been refined enough that it provides a suitable, albeit lower resolution, alternative to seismic monitoring for those periods when SAGA-BB is not available (Fig. 3b). This subspace detector cannot distinguish between explosions or regular LPs, cannot track changes in waveforms to nearly the same degree as a local network, and the magnitude threshold is relatively high, where in general any signal less than 8 m/s (as recorded at SAGA-BB) would not be detectable at regional distances. However, the primary advantage of the regional network is its more uniform homogeneity in time, which is important for not only operational monitoring in near real time but also for retrospective studies such as this one where the regional stations allowed to expand the seismic record of Sangay backwards in time, allowing to quantify the number of discrete seismic events since August 2007 (Fig. 3b). Interestingly, on a daily basis, regional detections could be higher than the quake counting in SAGA-BB. This is mainly due to the incomplete recordings from SAGA-BB (Fig. 2a). for waveform coherence relative to the original set of 5 basis functions across the 9 regional sensors. A subtle drop in the value of the sufficient statistic is seen towards the end of 2019 and into 2020, coincident with results obtained from the SAGA analysis, b) comparison between the number of seismic events from Sangay as detected at SAGA-BB (orange) and the nine regional seismic stations (blue). Red dots are the number of regional stations used to detect each discrete event. Green background shows quiescent phases (QP) while yellow, light orange and red backgrounds correspond to eruptive phases.

Table 1 summarizes the seismic activity at Sangay between August 2007 and May 2020 based on the regional network. Three main eruptive periods were identified. During the period from 4 August 2007 to 29 May 2013 the activity at Sangay was long-lasting and fluctuating. Average transient events reached 1817 events/d with a maximum of 72 and a median of 15 events/d. In contrast, between 2015 and 2018, four clear short episodic eruptions occurred one each year, as previously identified on SAGA-BB (Fig. 2b).

During this period, the four phases were statistically similar. They lasted ~100 days each and their transient events significantly increased in comparison with the previous period reaching ~8130 per day with a median of 86 and a maximum of 148 events/d. In addition, the 2019-2020 period showed fluctuating and long-lasting activity similar to what was observed during the first period between 2007 and 2013, but the daily event rate is three times higher. Average transient events reached 5250 events/day with a median of 40 and a maximum of 207 events/d, which corresponds to the maximum registered during the entire study period. Finally, the five repose periods recognized in between the eruptive phases display a minimum seismic activity with an average lower than one seismic event per day (Fig. 3). This agrees with the no escalating unrest before eruptions observed at SAGA-BB (Fig. 2b). Interestingly, the daily seismic event counting reconstructed for Sangay fits nicely with the satellite-based observations, especially with the thermal and VAAC alerts, validating and providing complementary data to understand eruptive dynamics (see section 3.3). This confirms the overall long-lasting activity since at least 2007 until 2013, the episodic activity from 2015 to 2018, and the increased eruptive activity starting in 2019 and still ongoing until the time of writing of this manuscript (March 2022).

SO 2 emissions

The ScanDOAS instrument worked intermittently. Nevertheless, it could record data from: the quiescent phase of 2013-2015 (52%), the 2015 eruption (72%), part of the quiescent phase between 2015 and 2016 (16%) and the ongoing eruption from November 2019 to May 2020 (23%). Percentages refer to the functioning days of the station in respect to the duration of the quiescence or eruptive phase.

During the first quiescence phase (2013)(2014)(2015), SAGA-DOAS detected between 1 and 4 valid measurements for 11-72 minutes of daily valid measurement duration (DVMD), which corresponds to the measuring time (in minutes) of all valid scans throughout a day (Fig. 4a, Table 1). The observed SO2 flux ranged from 2 to 219 t/d with an average of 1126 t/d (Table 1). During the 2015 eruptive period the station registered a slight increase in the number of valid measurements (1-6) and DVMD (18-108 minutes, Fig. 4a). However, the SO2 flux was 1032 t/d average with a maximum single measurement of 224 t/d, like the previous quiescence phase (Fig. 4b). In the second quiescence period (2015-2016): 1 to 9 valid measurements were detected with 10-162 minutes of DVMD. The SO2 flux reached 3953 t/d, with a maximum single measurement of 242 t/d (Table 1). Finally, the current eruption shows a significant increase in all the parameters. The daily valid measurement duration lies between 13-414 min with 1-28 valid measurements per day. The SO2 flux ranges from 12 to 2435 t/d with an average of 229277 t/d (Table 1). In summary, we identified an increase in DVMD and the number of valid measurements between quiescent and eruptive phases (Fig. 4a). Interestedly, the SO2 fluxes lay in the same order of magnitude during the 2013-2015 period, which corresponds to two quiescence (QP1 & QP2) and the 2015 eruptive phase (Table 1). In contrast, the 2019-2020 eruptive period displays a significant increase in SO2 measured degassing as well as daily valid measurement duration (Fig. 4b, Table 1). (Fig. 4b). ). During those eruptions, weekly volcanic radiative power was still lower than 60 MW (Fig. 5a). In contrast, the 2019-ongoing activity marked a significant increase in eruption intensity, reaching the highest weekly volcanic radiative power, up to 300 MW, with several minor peaks. The volcanic radiative energy curve also shows a significant increase in slope during this eruptive period (Fig. 5a). Figure 5b displays volcanic ash cloud occurrence according to Washington VAAC alerts. Plume heights are less than 3.5 km above the crater level, being similar throughout the entire period. Single isolated explosions in between eruptive phases also occurred. Throughout the study period, most of the ash clouds showed a westward drift (Fig. 5b), which agrees with the predominant wind direction in the Ecuadorian Andes. Finally, the cumulative number of ash alerts shows a significant increase during the 2018 and 2019ongoing eruptive phases (Fig. 5b).

Thermal power and ash plume height

Between 2001 and 2005 the activity was characterized by explosions and lava flows extruded from the Central Crater (Fig. 5c). Nonetheless, since 2007, ash emissions were focused on the central crater and the lava flows were emitted from the new Ñuñurcu vent (Fig. 5e &5g). Pyroclastic density currents related to lava front collapses also occurred throughout the entire study period, usually reaching a run-out distance of < 3 km towards the eastern and southeastern flanks (Fig. 5d &5f). However, during 2019-ongoing eruption, pyroclastic density currents reached more than 4 km from the vent, with a maximum length of 7.6 km on 11 December 2019 (Fig. 5h), which is significantly larger than those observed in any previous eruptive phase (Fig. 5d &5f). As a results, the occurrence of secondary lahars triggered by heavy rainfall that remobilizes this loose material has become a typical phenomenon during the ongoing eruption and has temporarily dammed the Upano river at its confluence with the Río Volcán, 25 km away from the summit (Fig. 1c).

Ash composition

Whole-rock analyses of bulk ash show a slight shift in silica content through time. This is plotted on the classification diagram of [START_REF] Peccerillo | Geochemistry of Eocene calc-alkaline volcanic rocks from the Kastamonu area, northern Turkey[END_REF]

Time-Averaged Discharge Rate (TADR) and cumulative volume estimation

The volumes of erupted lava and the TADRs were calculated by means of the thermal data acquired by satellite. The thermal approach is based on the empirical relationship between radiant energy (VRE) and erupted lava volumes (Vol) which can be expressed through a best-fit coefficient so that VRE/Vol = crad (radiant density) [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF]. It has been shown that this coefficient varies from volcano to volcano and is strongly influenced by the rheology of the erupted lava, the latter represented roughly by the silica content [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF]. We calculated the radiant density (crad) for Sangay using the silica content from the ash sampled at the SAGA-BB site and from literature data (Fig. 6). For the eruptive period between 2001 and 2013 we used the average 59.44wt% SiO2 for the 1990s eruptive products reported in [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF]. For the 2015, 2016 & 2017 eruptions we used 60.72wt% SiO2, which is the silica content of 2015 ashes, for the 2018 eruption 61.17wt% SiO2, and for 2019 and 2020 we used 59.6 and 59.0 wt% SiO2, respectively. The derived radiant densities are: 2.2910 7 (2001-2013), 1.8310 7 (2015-2017),

1.7010 7 (2018), 2.2210 7 (2019) and 2.4710 7 (2020) ±50% J/m 3 . Afterwards, we calculated minimum, mean and maximum time-average discharge rate and cumulative volume including uncertainties. It should be noted that the approach is valid under the fundamental assumption that radiant heat detected by MODIS is exclusively sourced from "effusive/extrusive activity" [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF]). However, if additional heat sources are produced by other volcanic phenomena/products (i.e. explosions, collapses of lava flow front, pyroclastic density currents) the volumes obtained with the thermal proxy are probably overestimated.

Since at Sangay volcano we are unable to quantify the proportion of heat sourced by "non-effusive activity" our estimate should be taken with caution and considered as maximum values. time-average discharge rate of 0.28±0.29 m 3 /s and a maximum of 2.07 m 3 /s (Table 1). During the episodic period a cumulative volume of 53.5±26.8 million m 3 was extruded with the mean TADR ranging from 0.02 to 3.77 m 3 /s. The 2015 eruption extruded 6.8±3.4 million m 3 , while the 2016, 2017 and 2018 eruptive phases emitted similar volumes of 15.8±7.9 million m 3 , 15.8±7.9 million m 3 and 15.1±7.6 million m 3 , respectively (Fig. 7b, Table 1). Finally, the 2019-2020 period shows a significant increase in both timeaverage discharge rate reaching 5±4.1 m 3 /s with a maximum of 21.8 m 3 /s and a cumulative volume of 172±85.8 million m 3 of lava (Fig. 7c, Table 1), making this last period the largest and most energetic during the past 20 years (Fig. 7a). Time-average discharge rate during the last period is five times higher than the previous ones. This significant increase could explain the high recurrence of lava front collapses and associated co-pyroclastic plumes (Fig. 5b), in addition to the intense erosion of Volcán ravine observed at the southeastern flank. Finally, the cumulative volume indicates that Sangay volcano has become more active over time, which implies a sustained increase of the global average discharge rate (Fig. 7a).

DISCUSSION

Interpretations of eruptive activity based on clustering of multiparametric monitoring data

Based on the simple temporal clustering of the seismic catalog, thermal and ash emission alerts (Figs. 3b, 5a & 5b), we identified well-defined quiescence episodes between seven eruptive phases from 2001 to 2020 (Fig. 8). We grouped these various phases into three main eruptive periods: mild long-lived activity from 2001 to 2013, episodic eruptions between 2015 and 2018, and continuous and very intense activity since 2019 (yellow, orange, and red backgrounds, respectively in Fig. 8). SiO2 compositions of Sangay's products in the 1990s ranged from 59.28 to 60.38 wt% [START_REF] Monzier | Sangay volcano, Ecuador: structural development, present activity and petrology[END_REF]. In 2015 and 2018 ash compositions were slightly more acidic (60.7 -61.1 wt%, respectively), while the current activity expelled ashes with lower SiO2 contents (59.6 and 59.0 wt% SiO2, Fig. 8a), similar to that of the 1990s. Overall, it seems that long-lasting activity matches with the less evolved magmas, while episodic activity takes place whenever compositions are slightly more acidic. The lower SiO2 content of the ongoing eruptive products coupled to a higher MgO content (2.63-2.87 wt% in Fig. 6d) indicates an injection of a more basic magma into the system. Given the few data and the quality of the ash samples, future work should focus on robust sampling of fresh lava to improve the geochemical patterns presented here and to better understand Sangay's magmatic system including magma reservoir(s) location, crystallization and magma replenishment processes. However, despite the low number of samples, the observed changes in SiO2 lie beyond the analytical error, indicating that they are meaningful.

SO2 degassing at Sangay started to be monitored at the end of 2013 (Fig. 8b green bars). For the period 2013-2015 daily SO2 fluxes are in the same order of magnitude (10 2 t/d) throughout both, quiescence and episodic eruptive phases. We interpret this as reflecting either, or a combination of: (1) crystallizationinduced volatile saturation [START_REF] Tait | Pressure, gas content and eruption periodicity of a shallow, crystallising magma chamber[END_REF]) of the remaining non-eruptible 2001-2013 magma, promoting its evolution towards more acidic compositions and its eruption by second boiling [START_REF] Putirka | Down the Crater: Where Magmas are Stored and Why They Erupt[END_REF]); (2) small magma injections, which only provided gas and heat that reactivated the shallow reservoir without reaching the surface [START_REF] Humphreys | Shallow-level decompression crystallisation and deep magma supply at Shiveluch Volcano[END_REF]). These two mechanisms (individually or combined) would cyclically have taken place triggering the four short-lived phases observed between 2015 and 2018, as well as the shift to slightly more acidic compositions, accompanied by similar SO2 emissions during quiescence and eruptive phases. Unfortunately, the station stopped working from 2015 to 2019, so there is not enough data to confirm if degassing levels stayed similar during quiescence and eruptive phases and to assess which mechanism(s) controlled the episodic eruptive activity during this period. Therefore, the observed episodic activity could also be explained by (3) cycles of crystallization-driven plug formation which lead to progressive gas accumulation triggering the episodic eruptions once the gas pressure overcame the resistance of the plug [START_REF] Battaglia | Autopsy of an eruptive phase of Tungurahua volcano (Ecuador) through coupling of seismo-acoustic and SO2 recordings with ash characteristics[END_REF][START_REF] Woitischek | The control of magma crystallinity on the fluctuations in gas composition at open vent basaltic volcanoes[END_REF]. Finally, for the 2019-ongoing eruption, both ground-based and satellite SO2 monitoring systems registered a significant increase in SO2 emissions (Fig. 8b). This could be interpreted as the refill of the shallow reservoir by a more basic and gas-rich magma that ascended from depth [START_REF] Putirka | Down the Crater: Where Magmas are Stored and Why They Erupt[END_REF]).

The regional seismic network allowed to identify low, oscillating, and long-lasting activity from 2007 to 2013 (Fig 8c). This fluctuating behavior could be related to gas waves that promoted episodes of more intense activity alternated with repose intervals over periods of days as proposed by [START_REF] Michaut | Eruption cyclicity at silicic volcanoes potentially caused by magmatic gas waves[END_REF].

On the contrary, between 2015 and 2018 intense but short-lived eruptions occurred, as can be observed in the increased number of explosion quakes per day (Fig 8c). This activity could have been triggered by degassing-induced second boiling during its respective preceding quiescence periods, as the SO2 data suggests. In 2019-2020 we observed an oscillating but increased number of quakes per day compared to the entire 2007-2018 period (Fig 8c). This could be related to a magma intrusion that ascended from depth, supporting the hypothesis presented previously based on the geochemical and SO2 data (Fig 8a & 8b).

Additionally, the occurrence of a single seismic family (A) from January 2015 to March 2020 suggests that there is only one dominating eruptive process (explosion type) throughout all active periods (Fig. 2). After

March 2020, two more seismic families appeared (C, D in Fig. 2), possibly indicating the involvement of new eruptive processes and/or significant changes in Sangay's morphology.

Eruptive plumes at Sangay are emitted via minor explosive activity and ash venting from the Central Crater.

The average height of the plumes has been of 1±0.6 km acl (Fig. 8d), corresponding to 6.3±0.6 km asl throughout the entire study period. Even though in 2019-2020 eruptive activity increased in various parameters, including the cumulative number of Washington VAAC alerts, the height of the eruptive plumes did not (Fig. 8). A possible explanation for the similar plume heights could be that explosions were triggered by slugs of gas that reached a certain threshold [START_REF] Woitischek | The control of magma crystallinity on the fluctuations in gas composition at open vent basaltic volcanoes[END_REF]). On the contrary, the high number of VAAC alerts is most likely associated to the increase of explosion frequency (given that there is more available gas), co-pyroclastic plumes triggered by various lava front collapses (Figs. 5d,5f & 5h) and/or the intense erosion (lateral, vertical, and down-and upstream) of Volcano ravine at the southeastern flank.

At Sangay, most of the thermal alerts are related to lava flow emplacement. In 2001-2013, the maximum time-average discharge rate (TADR) was 2.07 m 3 /s, in 2015 it reached 2.95 m 3 /s, in 2016: 5.65 m 3 /s, in 2017: 5.16 m 3 /s, in 2018: 3.55 m 3 /s and in 2019-2020 it increased drastically to 21.8 m 3 /s (Table 1). Overall, from 2001 to 2013 max TADR was the lowest throughout the entire study period (Fig. 8e), which could be related to the emptying of the shallow reservoir via minor but long-lasting (i.e. continuous) activity. In contrast, in the episodic period (2015)(2016)(2017)(2018) there was an increase in max time-average discharge rate from

2015 to 2016, while in 2017 it started to drop until 2018, forming a bell-shaped pattern (Fig. 5b & 8e). The 2001-2018 period could be interpreted as starting with a steady-state reservoir, which fed the system from 2001 to 2013. Then, between 2015 and 2018, second boiling and/or injections of new magma provided the necessary gas to trigger the episodic period that expelled the remaining magma of 2001-2013. Finally, the much higher max time-average discharge rate in 2019-2020 (Fig. 8e) is most likely related to a new voluminous and gas-rich injection of magma that ascended from depth, as was interpreted based on aforementioned monitoring parameters.

Similar changes in the eruptive behavior from long-lasting to episodic eruptions have been reported for the long-term activity of other volcanoes (e.g., Tungurahua, [START_REF] Hidalgo | SO2 degassing at Tungurahua volcano (Ecuador) between 2007 and 2013: Transition from continuous to episodic activity[END_REF]. Nevertheless, a detailed comparison with other similar systems is not possible given the lack of continuous data on Sangay, as well as the lack of long-term research studies published for open-system volcanoes.

Model of magma discharge based on injection, cooling and viscous relaxation time relationship

Given the long-lasting activity recorded pre-2013 and since 2019, we can argue that a persistently replenished shallow magmatic reservoir must be feeding the surface eruptive activity under a steady-state behavior, similar to that observed at basaltic open systems like Stromboli and Etna (Wadge and Guest 1981; [START_REF] Harris | Magma budgets and steady-state activity of Vulcano and Stromboli[END_REF]. The different intensity of the activity, which was mild pre-2013 and intense since 2019, could be controlled by the injection size and time. It seems that more voluminous basic and gas-rich magma intrusions with shorter injection times are responsible for the 2019-ongoing activity (Fig. 8), making this phase last longer and show more intense surface activity than before. Steady-state behavior like the one observed during these two periods could be recurrent at Sangay, nevertheless, we only have insights into the last 20 years of activity.

In between these two long-lasting eruptive periods, for 4 years, from 2015 to 2018, Sangay displayed episodic activity of shorter duration and milder surface phenomena. Silica content seems to slightly evolve from more basic magmas sampled in the 1990s to a slightly more differentiated composition for 2015 and 2018 (Fig. 8a). This could be in agreement with a progressive magma depletion and passive degassing of the shallow reservoir, followed by a crystallization induced heating and volatile exsolution of the upper remaining magma ending in eruption by second boiling [START_REF] Villareal R ; Morona | Plan estratégico institucional del Gobierno Autónomo Descentralizado del cantón Morona. Gobierno Autónomo Descentralizado[END_REF][START_REF] Tait | Pressure, gas content and eruption periodicity of a shallow, crystallising magma chamber[END_REF][START_REF] Putirka | Down the Crater: Where Magmas are Stored and Why They Erupt[END_REF]. Smaller intrusions could also have triggered the remaining slightly more differentiated magma to become eruptible again and produced the short duration episodes observed between 2015 and 2018.

We propose a simplified model based on [START_REF] Degruyter | A model for eruption frequency of upper crustal silicic magma chambers[END_REF] and [START_REF] Mittal | Volatile Degassing From Magma Chambers as a Control on Volcanic Eruptions[END_REF] to explain the geophysical and subaerial observations from Sangay. The model defines characteristic timescales for different processes and depending on the ratios between them, different regimes occur, either locking the magma (no eruption) or leading to an eruption (eruptions triggered by second boiling, mass injections or buoyancy) (Fig. 9). The timescales, in the more general model, are related to: magma injection, reservoir cooling, viscous relaxation of the crust, fluid pore pressure and its diffusion [START_REF] Degruyter | A model for eruption frequency of upper crustal silicic magma chambers[END_REF][START_REF] Mittal | Volatile Degassing From Magma Chambers as a Control on Volcanic Eruptions[END_REF]. These parameters are calculated from structural and rheological properties of the magmatic system (e.g., magma density and viscosity, size of reservoir, length of conduit, permeability of crust, thermal and pore diffusivities, etc.), most of which are unknown for Sangay.

However, the known surface activity, rest outgassing and/or erupting, can be used to locate any specific activity in Fig. 9 and propose a range of values for the ratio of timescales at this specific time. We assume that there were no significant changes over the last two decades regarding the geometry of the volcanic system as well as the viscous relaxation and cooling times. In this case, changes in the injection time, controlled by variations in magma inflow rate to the reservoir, result in a proportional change in the parameter space defined by relaxation to injection times vs. the cooling to injection times (Fig. 9).

Moreover, the fact that gas emissions occurred even during periods when other variables remained low, indicates that the reservoir is open to degassing, i.e., that the system maintains a non-zero permeability (efficient volatile diffusion), resulting in a model that combines the Degruyter and Huber (2014) diagram and the third regime defined by [START_REF] Mittal | Volatile Degassing From Magma Chambers as a Control on Volcanic Eruptions[END_REF] for the evolution of Sangay's eruptive activity (Fig. 9).

We hypothesize that the different eruptive periods observed at Sangay over the last years obey to different scaling between these timescales. For instance, if the geometry of the system remains unchanged, the rate of magma injection into the reservoir and/or changes in permeability would lead to transitions between the various regimes (Fig. 9). The fact that the height of the ash plumes is confined to a narrow range over the last 20 years suggests that the eruption style is similar throughout the different periods, while the frequency of the events and duration of each eruptive phase stay variable. The simplest way to explain these variations is to attribute the primary role in the dynamics to the rate of magma intrusion (i.e. magma injection time) [START_REF] Degruyter | A model for eruption frequency of upper crustal silicic magma chambers[END_REF][START_REF] Mittal | Volatile Degassing From Magma Chambers as a Control on Volcanic Eruptions[END_REF]. A decrease in the injection time due to an increase in the mass inflow rate could lead to a situation in which the time for pore diffusion is much higher than the time for magma intrusion. In both long-lasting periods, eruptions could have been triggered by mass injection or by buoyancy (controlled by volatile content), depending on how efficient viscous relaxation of the edifice was with respect to magma intrusion. Based on their different chemistry, volume and degassing, the 2001-2013 eruptions were most likely fueled by buoyancy, while the 2019-ongoing ones probably result from mass injection (Fig. 9). On the other hand, the episodic period between 2015 and 2018 could respond to a lower mass inflow rate (i.e. higher injection time), due to which the system produced sporadic eruptions by second boiling or remained quiet due to locking of the magma (Fig. 9). In all cases, the frequency of the eruptions is directly related to the mass inflow rate, which is a way to understand the changes from episodic to long-lasting activity at Sangay volcano.

CONCLUSION

Sangay is the most active volcano in continental Ecuador. For the last 20 years, we identified seven eruptive phases clustered into three major eruptive periods based on ground-based and satellite-derived information.

These periods were of long-lasting (2001-2013 & 2019-ongoing) and episodic nature (2015, 2016, 2017 & 2018). Surface activity remained similar over time; eruptive columns reached 1±0.6 km above the crater All the recorded parameters are consistent with open-vent activity explained by a variation in the deeper magma inflow rate into the reservoir, assuming that structural and rheological properties of the magmatic system remain confined to a relatively narrow range during the study period. Observation of these changes would not have been possible without the long-term multi-parametric monitoring presented here. These techniques, if adopted elsewhere, have the potential to provide new insights into the "unpredictable" behavior of other andesitic open system volcanoes worldwide.

  explored seismic and acoustic signals of the degassing explosions at Sangay recorded in April 1998, based on a temporary experiment conducted at a close range from the summit (<2 km). They distinguished three types of explosive events: (1) simple impulses (<1 min in duration) generated by degassing explosions at the vent, (2) extended degassing 'chugging' events (1-5 Hz and duration of 2-5 min) which could be related to resonances in the conduit, and (3) extended degassing events (>5 Hz and >5 min). Ortiz et al. (2020) showed near-constant activity from mid 2006 to mid 2012, with eight detections per day on average for Sangay based on a temporary infrasound array. Carn et al. (2008) reported negligible SO2 emissions at Sangay from 2004 to 2006 based on Total Ozone Mapping Spectrometer (TOMS) and Ozone Monitoring Instrument (OMI) satellite sensors. Morales Rivera et al. (2016) found consistent motion towards the satellite on the southwestern flank at 1.3±0.4 cm/y and away on the southeastern flank at 4.3±0.4 cm/y based on ALOS Interferometric Synthetic Aperture Radar (InSAR) ground deformation surveys between 2007 and 2011. Given the satellite's view direction, they suggested that the observed deformation was consistent with downhill motion of the steep flanks of the volcano by 3-5 cm/y.

Fig. 2

 2 Fig. 2 a) Functioning of SAGA-BB in percentage of time per day, b) Number of seismic events per day, black bars represent all detected events, while orange, red, green, and blue bars display the number of events belonging to each seismic family, c) amplitude of the events also separated by family and expressed in m/s. Green background shows quiescent phases (QP) while orange and red backgrounds correspond to eruptive phases. d) Similar waveforms showing monthly stacks for the main families. Each waveform corresponds to a stack of all similar events identified during a given date listed at the end of each line. The waveforms were detected during all periods of observation.

Fig. 3

 3 Fig. 3 a) Values of sufficient statistic of the subspace detector as a function of time, which can be considered a proxy

Fig. 4

 4 Fig. 4 SO2 data from the ScanDOAS station located at 6 km downwind from Sangay crater. a) Dark and light blue bars show the daily validated measurement duration (DVMD) and number of valid measurements respectively, b) comparison between daily SO2 fluxes (DOAS sensor) from December 2013 to May 2020 and the SO2 mass derived from (MOUNTS, TROPOMI sensor). Note that SO2 is expressed in t/d and masses derived from TROPOMI in tons. Gray backgrounds are periods without data from the DOAS instrument, green background shows quiescent phases (QP) while orange and red backgrounds represent eruptive phases. The satellite OMI sensor indicates negligible SO2 between 2004 and August 2019, while from September 2019 to May 2020, it reported sporadic SO2 emissions ranging from 30 to 1000 tons at Sangay. Additionally, TROPOMI data processed by MOUNTS identified daily SO2 degassing since the onset of the

Figure 5

 5 Figure 5 summarizes the observations of the subaerial activity of Sangay volcano between January 2001 and May 2020 recorded by various satellite instruments. Both thermal and ash emission alerts clearly show eruptive phases and quiescence over time (Fig. 5a & 5b). From 2001 to 2013, the eruptive activity was sparce, but long-lasting and the weekly Volcanic Radiative Power was lower than 60 MW (the maximum registered in this period). The activity paused between 2005 and 2006. After this activity interruption, in 2007, the Ñuñurcu vent appeared 210 m southeast of the Central Crater increasing the slope of the Volcanic Radiative Energy (Fig. 5a). Eruptive activity paused again in 2013 for a period of 21 months. Between 2015 and 2018, episodic eruptions occurred at the rate of 1/year, all of them separated by well-defined quiescent phases (Fig. 5a & 5b). During those eruptions, weekly volcanic radiative power was still lower than 60 MW

Fig. 5 a

 5 Fig. 5 a) Composite graph of the thermal anomalies recorded by MIROVA. On the left axis is the Volcanic Radiative Power (VRP) in megawatts and on the right axis weekly Volcanic Radiative Energy (VRE) in Joules. b) Washington VAAC confirmed alerts from January 2004 to May 2020 in meters above the crater level (m acl), cumulative number of alerts is shown on the right axis. Wind rose diagram displaying predominant wind direction and speed obtained from VAAC. Single explosions in between periods are shown as "S. Exp". Colored backgrounds separate the eruptive from quiescent phases (QP). Data downloaded from http://www.mirovaweb.it/ and

  in figure 6. Sangay's products are low-K basalts and andesites. Bulk tephra samples from 2015 to 2018 show similar ~61wt% SiO2 content, while 2019 and 2020 ashes are slightly more basic (59.6 -59.0wt% SiO2), closer to the 1990s magma compositions taken from Monzier et al. (1999) (Fig. 6a & 6b). MgO compositions also show an increasing trend confirming the similarity of the most recent products to the samples erupted during the 1990s (Fig. 6c & 6d).

Fig. 6

 6 Fig. 6 a) K2O vs. SiO2 diagram for Sangay's bulk ash compositions for 2015, 2018, 2019 and 2020 samples compared to previous geochemical data from Monzier et al. (1999). b) Zoom on K2O vs. SiO2 diagram showing the recent products compositions. c) MgO vs. SiO2 diagram, and d) Zoom on the MgO vs. SiO2 diagram. Error in analyses is smaller than the size of the symbols.

Fig. 7

 7 Fig. 7 TADR (black lines) and cumulative volume (red lines) for Sangay volcano derived from MIROVA thermal data and the methodology proposed by (Coppola et al. 2009, 2013, 2017). Colored backgrounds represent the eruptive stages and quiescence phases (QP) from a) January 2001 to May 2020, b) between 2015 and 2018, c) during 2019-2020 ongoing eruptive period. Note that black squares are the mean TADR and black and dashed red lines depict the error of the TADR and volume estimations, respectively.

Figure

  Figure 7a displays the estimated time-average discharge rate (TADR) and cumulative magma volume emitted throughout the last two decades. From 2001 to 2020, Sangay volcano extruded between 0.02 and 21.8 m 3 /s (TADR), which implies a cumulative extruded lava volume of 326±163 million m 3 . A third of the volume was emitted during the long-lasting period of 2001-2013, reaching 100±50 million m 3 with a

Fig. 8

 8 Fig. 8 Composite temporal series from 1990s to 2020. a) wt% SiO2 whole rock analysis of the eruptive products, b) daily SO2 flux (SAGA-DOAS) and mass (MOUNTS), c) quakes per day based on the regional seismic network, d) eruptive column height in km acl and cumulative Washington VAAC alerts, e) Time Average Discharge Rate (TADR) and cumulative volume of the extruded lavas based on VRP provided by MIROVA system. Colored backgrounds represent the eruptive (yellow, orange and red) and quiescence (green) phases. Dashed lines depict the onset of the different monitoring systems at Sangay.

Fig. 9

 9 Fig. 9 Simplified conceptual model of the various regimes that control the evolution of eruptive activity based on magma inflow rate into the reservoir and the relationship between injection, cooling and viscous relaxation times. Modified from Degruyter and Huber (2014) and Mittal and Richards (2019). Black dashed line indicates the possible trajectory in this model for Sangay eruptions throughout the last two decades, assuming negligible changes in the geometry of the volcanic system, the viscous relaxation and cooling times.

  level and lava flows of 1-3 km length were emplaced towards the eastern and southeastern flanks and frequently produced pyroclastic density currents due to lava front collapses. In contrast, magma composition has become slightly more basic over time, shifting from 60.7 wt% in 2015 to 59.0 wt% SiO2 in 2020, indicating new magma injections.Over the last 20 years, Sangay volcano has extruded at least 326±163 million m 3 of magma in three main periods. The first one emitted 100±50 million m 3 through long-lasting activity. The second expelled 54±27 million m 3 in four short-lasting eruptive episodes, which lasted ~100 days each. Finally, the 2019-ongoing eruptive period has been the most intense one in the last two decades. It has extruded more than 172±86 million m 3 (until 31 May 2020). The high discharge rate (TADR) in this period has triggered numerous lava front collapses (i.e., pyroclastic density currents) which have aggressively eroded Volcán ravine.Rainfall-triggered secondary lahars have at various times temporarily dammed the Upano river located 25 km downstream, marking the first time in the last two decades that lahars have been reported and have affected populated areas.

Table 1 .

 1 Statistical summary of the seismic and gas emission data based on 9 broadband regional seismic stations and the SAGA-DOAS station. Additionally, time-averaged discharge rate (TADR) and cumulative volume extruded at Sangay volcano are shown. *values are estimated from 2001. QP = Quiescent Phase, VM = valid measurements, DVMD = daily validated measurement duration.

	Eruptive	2007-	QP 1	2015	QP 2	2016	QP 3	2017	QP 4	2018	QP 5	2019-
	Phase	2013										2020
	Onset	04.08.07	30.05.13	04.01.15	19.04.15	08.03.16	31.07.16	20.07.17	28.07.17	11.08.18	09.12.18	06.05.19
	End	29.05.13	03.01.15	18.04.15	07.03.16	30.07.16	19.07.17	27.10.17	10.08.18	08.12.18	05.05.19	31.05.20
	Duration (d)	2125	584	105	324	145	354	100	287	120	148	392
	Average seism/d	18±17	0.1±0.3	72±30	0.3±0.8	88±25	0.5±0.9	83±39	0.1±0.4	80±28	0.6±1	52±50
	Median seism/d	15	0	77	0	95	0	93	0	86	0	40
	Max. seism/d	72	5	133	6	133	5	148	3	126	7	207
	Cum. Seism	38781	31	7568	100	12767	187	8330	40	9595	86	20339
	#VM SO2	-	1-4	1-6	1-9	-	-	-	-	-	-	1-28
	DVMD (min)	-	11-72	18-108	10-162	-	-	-	-	-	-	13-414
	SO2 flux (t/d)	-	11±26	10±32	39±53	-	-	-	-	-	-	229±277
	Max SO2 flux (t/d)	-	219	224	242	-	-	-	-	-	-	2435
	TADR (m 3 /s)	0.28±0.29*	-	0.8±0.56	-	1.3±1.16	-	1.35±1.08	-	1.47±0.63	-	5±4.1
	TADR max (m 3 /s)	2.07*	-	2.95	-	5.65	-	5.16	-	3.55	-	21.8
	Cum. Volume		-		-		-		-		-
	(Mm 3 )	100±50*		6.8±3.4		15.8±7.9		15.8±7.9		15.1±7.6		172±85.8

Acknowledgement

The first author thanks the late Luis Peña and Wilson Franco for their friendship and unconditional support to Vasconez family, which allowed them to begin and successfully end his educational career. The authors would like to thank the entire staff of ECU-911 Macas, in particular to Luis Castillo for gently sharing pictures and videos of Sangay eruptions. Special acknowledgment goes also to Byron Rivadeneira, Diego Barba, Jorge Anhalzer, Eric Brosier and Grupo Aéreo del Ejercito Ecuatoriano for their cooperation in station maintenance. We also thank the entire staff of Instituto Geofísico of Escuela Politécnica Nacional (IG-EPN) in particular to Marco Antonio Solis, Jéssica Mejía, Paola Nazate, Luis Velez, Daniel Cárdenas, Sandro Vaca, Freddy Vásconez, Francisco Mejía, Iván Tapa, Ana Peralvo and Andrea Cordova. We also thank Pablo Samaniego, who gently shared the geochemistry data from Monzier et al., 1999. Finally, we thank the two anonymous reviewers and the editor for their insightful and constructive comments which helped to improve the original manuscript. Part of this investigation was funded thanks to the PAPIIT project IA102221. This research was conducted in the context of IG-EPN's project "Generación de Capacidades para la Emisión de Alertas Tempranas" funded by Secretaría Nacional de Planificación y Desarrollo (SENPLADES) and the Laboratoire Mixte International Seismes et Volcans dans les Andes du Nord (LMI SVAN).

provided the MIROVA data and helped in their analysis. SV provided de SO2 data and helped in their analysis. PR identified the change from

Data downloaded from http://www.mirovaweb.it/ and https://www.ssd.noaa.gov/VAAC/messages.html respectively. From c) to h) composition of photographs and thermal images illustrating the active vents and products over the study period and maps displaying the locations of MODIS and VIIRS thermal anomalies centers from the different eruptive phases based on FIRMS https://firms.modaps.eosdis.nasa.gov/map/.

Data availability

Supporting information is available alongside the online version of this manuscript. Supplementary Material 1: Subspace detector and binary classifier and Supplementary Material 2: Geochemical data.


Raphaël Chattot
email: raphael.chattot@umontpellier.fr

Marta Mirolo

Isaac Martens

Kavita Kumar

Vincent Martin

Amir Gasmi

Laetitia Dubau

Frédéric Maillard

Luis Castanheira

Jakub Drnec

Beware of cyclic voltammetry! Measurement artefact in accelerated stress test of fuel cell cathode revealed by operando X-ray diffraction

Keywords: Nanocatalysts, Fuel Cells, Accelerated Stress Test, Cyclic Voltammetry, Operando X-ray Diffraction

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

TEXT 1. Introduction

Cost and durability remain two major barriers to the widespread deployment of proton exchange membrane fuel cell (PEMFC) systems, especially regarding fuel cell electric vehicle application. Because the degradation of the costly platinum group metal(s) (PGM)-based cathode catalyst accounts for most of the system cost and performance decay, a lot of efforts are dedicated to design more efficient and durable cathode catalyst layers [START_REF] Kodama | Challenges in Applying Highly Active Pt-Based Nanostructured Catalysts for Oxygen Reduction Reactions to Fuel Cell Vehicles[END_REF] . To be competitive, practical PEMFC cathode should (among other prerequisites) exhibit >0.44 A mgPGM -1 at 0.9 V cell voltage, with at least 60 % activity retention over 8,000 h of operation (2025 targets from the Fuel Cell Technical Team of the U.S Department of Energy, DoE/FCTT) [START_REF] Abdel-Baset | The Fuel Cell Technical Team Roadmap[END_REF] . To shortcut the duration of such long durability testing by at least two orders of magnitude, accelerated stress tests (ASTs) have been designed to investigate the viability of commercial and experimental catalyst materials in membrane electrode assemblies (MEAs). ASTs typically consist in applying potential/current cycling profiles to the MEA, simulating in a high frequency fashion the electrochemical stress encountered during operation. Despite numerous AST protocols can be found in the literature (including different potential profiles, potential sweep rate, different upper and lower potential limits etc. [START_REF] Stariha | Recent Advances in Catalyst Accelerated Stress Tests for Polymer Electrolyte Membrane Fuel Cells[END_REF][START_REF] Zhang | The Impact of Potential Cycling on PEMFC Durability[END_REF][START_REF] Borup | Recent Developments in Catalyst-Related PEM Fuel Cell Durability[END_REF][START_REF] Uchimura | Electrocatalyst Durability under Simulated Automotive Drive Cycles[END_REF] ) the DoE/FCTT proposed 30,000 square-wave cycles made of potential steps between 0.6 V (3s) and 0.95 V (3s) as a standard protocol (MEA notably operated at 80 °C with hydrogen and nitrogen supplied at the anode and cathode, respectively). However, for the impact of such AST on the cathode catalyst electrocatalytic performance (and to some extents structural) stability to be possibly evaluated, the protocol also includes beginning of test (BoT), end of test (EoT) and intermediate electrochemical characterizations. Such characterizations consist in cell voltage and electric current measurements under reactive gas (polarization curves for activity estimation) and cyclic voltammetry (CV) experiments under inert atmosphere.

Cyclic voltammetry is a powerful tool used as routine technique in numerous fields of electrochemistry, as it probes the redox events occurring at the electrode-electrolyte interface. In case of Pt-based electrodes, measuring the desorption coulometry of under-potentially deposited hydrogen (HUPD) allows (supposing the maximum surface density of hydrogen adsorbed a known constant) to evaluate the electrochemically active surface area (ECSA) of the electrode. This only requires short time potential excursion(s) of the electrode close to 0 V vs. the reversible hydrogen electrode (RHE), typically 0.05 or 0.1 V vs. RHE. Because the Pt-based PEMFC anode counter electrode in hydrogen atmosphere can be considered as an RHE under simple experimental conditions, CV experiments can be easily performed directly in the fuel cell device.

Supported by the stability of Pt metal observed in the Pourbaix diagram potential-pH range [START_REF] Pourbaix | Atlas of Electrochemical Equilibria in Aqueous Solutions[END_REF] , this method is considered as non-destructive. Consequently, CV experiments are an essential part of any AST protocol, since to pass the stability test, the catalyst should retain at least 60% of its initial activity and ECSA according to the DoE/FCTT. This implies ECSA has to be measured by CV experiments at least at BoT and EoT, but additional measurements also after 10, 100, 1k, 3k, 10k and 20k cycles are also recommended.

However, if pure Pt metal is likely the state of the catalyst at the BoT, Pt oxidation and dissolution are wellidentified fundamental pillars of catalyst degradation during fuel cell operation and AST [START_REF] Stariha | Recent Advances in Catalyst Accelerated Stress Tests for Polymer Electrolyte Membrane Fuel Cells[END_REF][START_REF] Uchimura | Electrocatalyst Durability under Simulated Automotive Drive Cycles[END_REF][START_REF] Topalov | Towards a Comprehensive Understanding of Platinum Dissolution in Acidic Media[END_REF][START_REF] Fuchs | Structure Dependency of the Atomic-Scale Mechanisms of Platinum Electro-Oxidation and Dissolution[END_REF][START_REF] Lopes | Eliminating Dissolution of Platinum-Based Electrocatalysts at the Atomic Scale[END_REF][START_REF] Sandbeck | Particle Size Effect on Platinum Dissolution: Practical Considerations for Fuel Cells[END_REF][START_REF] Kasian | Electrochemical On-Line ICP-MS in Electrocatalysis Research[END_REF][START_REF] Ehelebe | Evaluating Electrocatalysts at Relevant Currents in a Half-Cell: The Impact of Pt Loading on Oxygen Reduction Reaction[END_REF] . Pt cations can spread within the catalytic layer, being possibly further reduced in situ onto metallic nanoparticles according to the Oswald ripening process, or in the polymeric membrane by the crossover hydrogen [START_REF] Borup | Scientific Aspects of Polymer Electrolyte Fuel Cell Durability and Degradation[END_REF][START_REF] De Bruijn | Review: Durability and Degradation Issues of PEM Fuel Cell Components[END_REF][START_REF] Dubau | A Review of PEM Fuel Cell Durability: Materials Degradation, Local Heterogeneities of Aging and Possible Mitigation Strategies[END_REF] .

Thus, considering the aged catalyst state being a combination of Pt in metallic, oxidized and ionic states, the question whether a CV measurement is still a non-destructive tool in regards of these new moieties may deserve more attention. In fact, the same Pourbaix diagram of Pt which supports the stability of Pt metal in a typical CV experiment potential range clearly predicts the reduction of both Pt cations and oxides at potential approaching 0 vs. RHE. Most importantly, such low potential excursion of the cathode is not supposed to occur in practical fuel cell operation.

In this contribution, we use a combination of operando X-ray diffraction and ex situ characterization techniques to investigate the degradation of a Pt/C cathode catalyst in X-ray transparent PEMFC device during a typical AST procedure. The results reveal that, despite the apparently consistent agreement obtained from all the characterization techniques, the degradations observed in typical AST experiment analysis approach are likely misinterpreted. The 'measurement artefact' in AST protocols, where the structure of the aged catalyst is in fact irreversibly modified by the auxiliary CV characterisation tool, is presented and discussed in light of additional fuel cell tests, on-line inductively coupled plasma mass spectrometry (ICP-MS) and in situ WAXS measurements in liquid electrolyte. A commercial 40 wt.% Pt/C catalyst (40 wt. % Pt on medium surface carbon black, ID-n° 5097373, lot EC200840) from Heraeus (Hanau, Germany) was incorporated into 5 cm 2 MEAs (see experimental details in the Supplementary Material). For operando X-ray diffraction experiments, the MEA was mounted inside an X-ray transparent PEM single cell, described in a dedicated contribution [START_REF] Martens | X-Ray Transparent Proton-Exchange Membrane Fuel Cell Design for in Situ Wide and Small Angle Scattering Tomography[END_REF] current density steps between 0 and 1 A cm -2 was performed, and the cell voltage was measured after 3 min of stabilization at each step. During the polarization curve, the O2 stoichiometry was maintained at 9.5, while the H2 flow was kept constant (104 sccm, corresponding to a stoichiometry of 3 at 1 A cm -2 ).

Experimental

Afterwards, O2 gas at cathode was replaced by 250 sccm N2, and BoT CVs (3 cycles between 0.05 and 1.0 V vs. RHE at a potential sweep rate of 50 mV s -1 ) were recorded. The core of the AST consisted in 10,000 square-wave cycles made of potential steps between 0.6 (3s) and 1.0 V vs. RHE (3s) under N2 cathode gas atmosphere. Importantly, for each WAXS measurement (frequency of the measurements being shown later)

the potential cycling was interrupted and the cathode was maintained at E = 0.4 V vs. RHE during WAXS acquisition (for ca 5 min). Finally, EoT characterizations consisted in EoT CV in N2 cathode atmosphere (as CV BoT) followed by EoT polarization under O2 cathode atmosphere (as BoT polarisation curve). The cell was eventually turned off, and the aged MEA kept for further post mortem analysis. The same experiment protocol was also conducted in a 'conventional' single PEMFC (see picture in Figure S1) at the laboratory scale without operando WAXS measurements (see the Supplementary Material). Finally, complementary on-line inductively coupled plasma mass spectrometry (ICP-MS) and in situ WAXS measurements in liquid electrolyte were also performed, as introduced later and described in the Supplementary Material). Because the X-ray transparent fuel cell is not optimized for high performance,

Results and Discussion

j (A cm -2 MEA )
Begining Interestingly, however, a theoretical Pt specific surface area loss of 56 % is predicted from the nanoparticle coherent domain size increase, in excellent agreement with the measured ECSA loss of 59 %. This suggests that 95 % of the ECSA loss can be accounted for particle size increase, with no or only few losses of Pt metal mass. These results are consistent with numerous previous reports on Pt/C degradation in PEMFC 14- 16,20 , which highlighted the limited time existence of Pt z+ ions in operating PEMFC MEAs due to various redeposition processes [START_REF] Dubau | A Review of PEM Fuel Cell Durability: Materials Degradation, Local Heterogeneities of Aging and Possible Mitigation Strategies[END_REF] . Regarding the amplitude in particle size increase, reaching ca 10 nm after degradation tests in PEMFC is also quite typical. Martens et al., [START_REF] Martens | Imaging Heterogeneous Electrocatalyst Stability and Decoupling Degradation Mechanisms in Operating Hydrogen Fuel Cells[END_REF] reached 10.3 nm from commercial Pt/C MEA in similar conditions after 10,000 triangular wave AST between 0.6 V vs. RHE and 1.0 V vs. RHE at 50 mV s -1 . Similar particle size values were obtained also with different experimental conditions, such as potential cycling to different upper potential limits, different amount of cycles or even after extensive operation at constant current load [START_REF] Borup | Scientific Aspects of Polymer Electrolyte Fuel Cell Durability and Degradation[END_REF][START_REF] De Bruijn | Review: Durability and Degradation Issues of PEM Fuel Cell Components[END_REF][START_REF] Dubau | A Review of PEM Fuel Cell Durability: Materials Degradation, Local Heterogeneities of Aging and Possible Mitigation Strategies[END_REF][START_REF] Weiß | Impact of Intermittent Operation on Lifetime and Performance of a PEM Water Electrolyzer[END_REF][START_REF] Shokhen | Impact of Accelerated Stress Tests on the Cathodic Catalytic Layer in a Proton Exchange Membrane (PEM) Fuel Cell Studied by Identical Location Scanning Electron Microscopy[END_REF] . We now present the WAXS data measured operando over the course of the complete AST protocol. WAXS patterns were collected in grazing incidence configuration as represented in Figure 1.c. Importantly, the full MEA was scanned vertically by the X-ray beam from the anode to the cathode, and the WAXS patterns from the cathode were averaged for a complete capture of the degradation over the catalyst layer thickness 56 % surface loss [START_REF] Schröder | Tracking the Catalyst Layer Depth-Dependent Electrochemical Degradation of a Bimodal Pt/C Fuel Cell Catalyst: A Combined Operando Small-and Wide-Angle X-Ray Scattering Study[END_REF] . WAXS was measured between each experimental step, and after different amounts of cycles during the AST. Figure 4 shows the evolution of the coherent domain size over the course of the complete AST protocol. Strikingly, if the coherent domain size values from operando WAXS measurements agree with the values found ex situ at BoT and EoT, they reveal the observed coherent domain size increase mainly originates from the cyclic voltammetry experiment performed after the AST (from 4.3 nm to 10.4 nm), and not from the AST itself (from 3.8 nm to 4.3 nm). Closer look on the AST cycles in Figure 4.b shows that the AST cycles did induce a sharp coherent domain size increase during the first 1,000 cycles, but was followed by a much slower evolution for the rest of the test (as also previously observed [START_REF] Martens | Imaging Heterogeneous Electrocatalyst Stability and Decoupling Degradation Mechanisms in Operating Hydrogen Fuel Cells[END_REF][START_REF] Martens | Decoupling Catalyst Aggregation, Ripening, and Coalescence Processes inside Operating Fuel Cells[END_REF][START_REF] Hornberger | Impact of Carbon and Pyridinic-N Content on the Fuel Cell Performance and Durability of Carbon-Supported Pt Nanoparticle Catalysts[END_REF] ). To discard any experimental artefact possibly altering the results observed during the operando WAXS experiment (electrode potential control loss or beam damage), we performed ex situ X-ray diffraction (XRD) from laboratory source (see details in the Supplementary Material) on pristine and aged GDEs in the conventional laboratory cell with and without EoT CV experiment after the AST. As shown in Figure 4.a, the Scherrer formula applied on these independent measurements (see Figure S3) reveals an evolution of the coherent domain size from 4.01 nm to 5.87 nm or 9.04 nm after AST or after AST and EoT CV experiment, respectively, in agreement with the operando WAXS. This unambiguously reveals the detrimental, and largely overlooked impact a simple CV experiment can have on an aged catalyst sample (note the CV recorded at the BoT had no such impact). We argue that such dramatic changes must be caused by the strong reducing conditions experienced by the aged catalyst layer once its potential is swept in the HUPD region for ECSA measurement. The critical potential limit triggering this process is unfortunately not known at this stage, but it must be below E = 0.4 V vs. RHE, considering the WAXS measurements in The first two dissolution peaks, referred as A1 and C1 correspond to well-known anodic and cathodic processes, respectively [START_REF] Kasian | Electrochemical On-Line ICP-MS in Electrocatalysis Research[END_REF] . A1 is triggered by the formation of Pt oxides when the potential exceeds ca 0.9

V vs. RHE and C1 by the reduction of Pt oxides when the potential is decreased below ca 1.0 V vs. RHE.

Interestingly, the present on-line ICP-MS data (but in fact also others [START_REF] Sandbeck | Particle Size Effect on Platinum Dissolution: Practical Considerations for Fuel Cells[END_REF][START_REF] Kasian | Electrochemical On-Line ICP-MS in Electrocatalysis Research[END_REF] ) clearly allows the identification of a second cathodic dissolution process, labelled C2, when the potential is further decreased below ca 0.35

V vs. RHE. By extension of the mechanism ruling C1, we argue the presence of C2 is an indirect detection of the reduction of Pt surface oxides occurring a such low electrode potential. Since C2 is about 40 times less intense than C1, the amount of such remaining surface oxides must be significantly low.

Synchronization of the on-line ICP-MS data with in situ WAXS data measured in similar conditions (but with a potential sweep rate of 20 mV s -1 ) displayed in Figure 5.c reveals C2 occurs just before and during the adsorption of hydrogen (monitored by the lattice strain increase [START_REF] Chattot | Electrochemical Strain Dynamics in Noble Metal Nanocatalysts[END_REF] ). Interestingly, the C2 process seems to be magnified during the first of the three cycles recorded, suggesting it is highly sensitive to the history of the electrode. In Figure 5, the electrode was equilibrated at E = 0.4 V vs. RHE before performing the measurements. In the case of post-AST, aged catalyst, the amount of Pt oxides is likely even higher. In fact, and as shown in Figure S4, a mild additional reducing current contribution can be observed on the cyclic voltammograms during the first excursion at low electrode potential only at EoT. The present hypothesis is also supported by a previous contribution from Sandbeck et al. [START_REF] Sandbeck | Particle Size Effect on Platinum Dissolution: Considerations for Accelerated Stability Testing of Fuel Cell Catalysts[END_REF] , which identified from online ICP-MS dissolution profiles the lower potential limit of 0.6 V vs. RHE being too much positive in typical AST procedure to fully reduce Pt oxides on 4 nm nanoparticles. If these extra measurements in liquid electrolyte support the possibility of highly oxophilic Pt sites reduction occurring only at low electrode potential, this can hardly explain on its own such particle size increase. In fact, based on a cuboctahedral shape model [START_REF] Montejano-Carrizales | Direct Enumeration of the Geometrical Characteristics of Clusters[END_REF] , Pt nanocrystallites of 4.3 nm and 10.4 nm (before and after EoT CVs) are composed of ca 3,000 and 40,000 atoms, respectively. This means, ca 92 % of the total Pt mass in the form of Pt z+ cations or amorphous oxide surface layer would be needed to produce such a nanoparticle coherent domain size increase after reduction/redeposition, which seems unrealistic. The last allowed to clearly identify the aggregation process as anterior to the coalescence process in the overall degradation mechanism. Moreover, hydrogen adsorption is well-known to enhance Pt surface reconstruction [START_REF] Horch | Enhancement of Surface Self-Disfusion of Platinum Atoms by Adsorbed Hydrogen[END_REF] and nanoparticle sintering under thermal treatment [START_REF] Fiedorow | The Sintering of Supported Metal Catalysts. II. Comparison of Sintering Rates of Supported Pt, Ir, and Rh Catalysts in Hydrogen and Oxygen[END_REF][START_REF] Dubau | Defects Do Catalysis: CO Monolayer Oxidation and Oxygen Reduction Reaction on Hollow PtNi/C Nanoparticles[END_REF] or in electrochemical environment [START_REF] Mccrum | Hydrogen-Induced Step-Edge Roughening of Platinum Electrode Surfaces[END_REF] . Thus, from the available literature and the data of this contribution, we propose the following simple mechanism to covertly 'corrupt' standardized AST protocols (Figure 6

CONCLUSIONS

In this contribution, we used operando X-ray diffraction to capture the dramatic impact of cyclic voltammetry experiments on the structure of a Pt/C catalyst when performed on degraded (aged) fuel cell catalyst layer. Our observations, confirmed by a replica of the experiment, reveal that the cyclic voltammetry measurement performed at the end of a typical AST procedure causes more coherent domain size increase than the AST itself. Thus, despite being a very convenient tool to probe the catalyst structure in situ, the admitted non-destructive property of this technique on pristine (metallic) catalysts does not seem to hold upon catalyst degradation. The 2-3-fold increase in coherent domain size observed after cycling to

Figure 1 :

 1 Figure 1: Experimental setup and protocol. a-c) X-ray transparent PEMFC for operando WAXS; d) timeline of the different experimental steps and characterization tools used during the catalyst lifetime.

shown in Figure 1 b

 1 and Figure S1. This cell allows the acquisition of wide-angle X-ray scattering (WAXS) patterns of the isolated cathode operando, in grazing-incidence configuration (Figure 1.c). The various experimental procedures experienced by the catalyst powder are summarized in the timeline schematic of Figure 1.d. After the Pt/C powder catalyst being integrated into a catalyst layer and mounted in the PEMFC, the cell was heated at T = 80 °C and 100 % relative humidity (R.H.), and the anode and cathode sides were fed with 104 sccm H2 and 250 sccm O2 gas, respectively, with 0.5 bar backpressure (1.5 bar absolute). The cell was conditioned via a short protocol (considering limited allocated beamtime) consisting in 5 repetitions of a potential cycle composed of 2 min at open circuit voltage (OCP), followed by 3 min at cell voltage of 0.85 V and 10 min at 0.65 V. Then, a BoT polarization curve consisting of 10

Figure 2 :

 2 Figure 2: BoT and EoT electrochemical data in the conventional laboratory fuel cell. a) geometric surface area normalized current-voltage polarization curves, b) Pt surface area normalized current-voltage polarization curves and c) cyclic voltammograms recorded at 50 mV s -1 . In c), the Pt surface area (ECSA) is calculated by HUPD coulometry. Cathode loading: 0.2 mgPt cm -2 .

Figure 2 .

 2 a-b shows the cell performance measured at BoT and EoT in the conventional cell (results from the X-ray transparent cell are displayed in Figure S2).

Figure 2 .

 2 a clearly reveals the detrimental impact of the AST on the catalyst -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1

Figure 3 :

 3 Figure 3: BoT and EoT ex situ characterizations. a) WAXS patterns with corresponding coherent domain size extracted from Rietveld refinement and Pt specific surface areas. TEM images with corresponding number-and volume-averaged particle size distributions at b) BoT and c) EoT, respectively. On the TEM images, top left corners in b) and c), circles scale with the volume-averaged distribution mean are displayed.

Figure 4 .

 4 Figure 4.b were performed at this potential, and as confirmed below.

Figure 4 :

 4 Figure 4: Coherent domain size extracted from Rietveld measurement of the WAXS patterns measured ex situ (black markers) and operando (red markers), volume-averaged particle size estimated by TEM ex situ (blue markers) and coherent domain size estimated with the Scherrer formula from ex situ laboratory XRD performed on independent PEMFC tests. Operando WAXS was measured between a) each experimental step and b) after different amounts of cycles during the AST.

Figure 5 :

 5 Figure 5: Online inductively-coupled plasma mass spectrometry (ICP-MS) and in situ WAXS data measured during cyclic voltammetry experiments in 0.1 M HClO4 electrolyte. The potential profile in a) is valid for both the on-line ICP-MS and WAXS data in b) and c), respectively. The red boxes are displayed to guide the eyes in the identification of the C2 dissolution process.

  above-mentioned causes of coherent domain size increase is the sudden coalescence of aggregates. In contrast to the other possibilities discussed, such scenario could easily explain the ca 2.3fold size increase observed here. In a recent contribution[START_REF] Martens | Decoupling Catalyst Aggregation, Ripening, and Coalescence Processes inside Operating Fuel Cells[END_REF] , Martens et al. monitored operando the particle size increase of a commercial Pt/C catalyst over the course of an AST in PEMFC with wide-and smallangle X-ray scattering (WAXS vs. SAXS). The complementary sensitivities of the techniques regarding the particle size definition (coherent domain size vs. geometric size for WAXS and SAXS, respectively)

  ): (i) the core of the AST (cycling between 0.6 ≤ E ≤ 1.0 V vs. RHE) indeed triggers the well-known Pt degradation processes (Pt dissolution, nanoparticles aggregation and growth according to the Ostwald ripening process, Figure 6.a), (ii) the surface state reached after the potential cycling above 0.6 V vs. RHE is still partially oxidized; (iii) electrode excursion to reducing potentials of the 'HUPD region' during ECSA measurement at EoT reduces Pt surface oxides while hydrogen adsorption enhance Pt atoms surface diffusion and coalescence of aggregates (Figure 6.b).Finally, if the present study reports a quite consequent structural transformation caused by the cyclic voltammetry experiment performed at the EoT, the coherent domain size increase must be correctly interpreted. The nanoparticles do not simply grow from ca 4 nm to ca 10 nm, but aggregates formed during the AST coalesce during the EoT CV experiment. If the impact of such coalescence on the ECSA value cannot be easily accessed (at least not by conventional CV experiment), polarization curves measured after the AST but before and after the EoT CV experiment in FigureS5reveal a possible minor positive impact of the coalescence on the cell performance in this case. Note, however, the conclusions could be different if the CV experiments are performed regularly after various amounts of cycles as recommended by the DoE/FCTT. In fact, Castanheira et al.[START_REF] Castanheira | Accelerated Stress Tests of Pt/HSAC Electrocatalysts: An Identical-Location Transmission Electron Microscopy Study on the Influence of Intermediate Characterizations[END_REF] reported that implementing intermediate characterizations (CVs and COads stripping experiments) during a 96 h potentiostatic hold at E =1.0 V vs. RHE AST results in a 3.5-fold magnification of ECSA loss compared to the AST alone.

Figure 6 :

 6 Figure 6: Proposed mechanism highlighting the structural and chemical changes occurring during a) the core of the AST and b) the EoT CV experiment. The only few amounts of Pt oxides represented after AST refers to the oxides still present at E = 0.4 V vs. RHE, and further reduced during the EoT CV experiment.

  suggests a lightning coalescence of aggregated nanoparticles, triggered by the reduction of Pt oxides. Since the cyclic voltammetry experiments in the context of ASTs are used to determine the electrochemically active surface area of the catalyst, the question whether this measurement remains valid considering the structural changes it causes is now open, and alternative catalyst performance/structural assessment metrics may be needed.

Table 1 .

 1 activity, as a negative shift of ca 35 mV at 1 A cm MEA -2 is observed in the polarization curve. The CVs recorded at BoT and EoT in Figure2.c, likely provide the cause of such performance decrease. From HUPD coulometry, the AST caused nearly 60% of ECSA loss. Interestingly, normalizing the cell performance by the ECSA in Figure2.b highlights an increase of the Pt specific activity after the AST. Such opposite trends in mass-(or geometric) activity vs. specific activity can be explained at least partially by an increase of the nanoparticle size during ageing[START_REF] Perez-Alonso | The Effect of Size on the Oxygen Electroreduction Activity of Mass-Selected Platinum Nanoparticles[END_REF] . Importantly, despite the performance being different from one cell to the other, an excellent agreement is obtained for the ECSA trend, as shown in Moreover, the Pt specific

			a)			of Test End of Test		Voltage (V)	0.8 1.0	b)	End of Test Begining of Test
			H 2 /O 2						Cell
			100% RH 80 °C 0.5/0.5 bar backpressure 9.5x stoech. O 2 Constant H 2 flow (3.0x stoech. @ 1 A cm -2 MEA )			iR-free	0.6
										0.4
												0	10	20	30	40
													j (mA cm -2 Pt )
		100	c)
		50			48.8
	(mA)	0
	i	-50
							59	%	ECSA	loss
					19.9
		-100											Begining of Test
													End of Test
		-150
			0.0	0.1	0.2	0.3	0.4		0.5			0.6	0.7	0.8	0.9	1.0	1.1
							iR-free Cell Voltage (V)

-1

measured during a control experiment in liquid electrolyte for the same catalyst (not shown).

Ensuring a near-complete electrode-electrolyte contact from the catalyst layer is essential for further consistent interpretation of such AST.

Table 1 :

 1 Evolution of the ECSA measured by HUPD coulometry at BoT and EoT in the two PEMFC used.More insights on the origin of ECSA loss can be obtained from preliminary and post mortem analyses of the materials aged in the X-ray transparent cell. Both the fresh catalyst powder and aged catalyst powder collected from the tested MEA were investigated with TEM and ex situ WAXS (see experimental details very close to the volume-averaged particle size found by TEM (from 3.8 nm to 13.0 nm). Note that only the volume-averaged TEM sizes can be tentatively compared to the volume-sensitive WAXS, and that coherent domain size is a slightly distinct notion from nanoparticle geometric size. Furthermore, estimating nanoparticle size on the aged catalyst accurately from very local TEM images is a difficult if not impossible task considering the presence of aggregates and broad size distribution. Supposing cuboctahedral shape for the crystallites, it is possible to estimate the Pt specific surface areas at BoT and EoT from the WAXS coherent domain sizes[START_REF] Montejano-Carrizales | Direct Enumeration of the Geometrical Characteristics of Clusters[END_REF] . The values of 66 m 2 𝑔 𝑃𝑡 -1 and 29 m 2 𝑔 𝑃𝑡 -1 for BoT and EoT, respectively, are significantly higher than the values of 47.2 m 2 𝑔 𝑃𝑡 -1 and 19.2 m 2 𝑔 𝑃𝑡,𝑖𝑛𝑖

		X-ray transparent PEMFC	Conventional PEMFC
	ECSA BoT	47.2 m 2 g 𝑃𝑡 -1	48.8 m 2 g 𝑃𝑡 -1
	ECSA EoT	19.2 m 2 g 𝑃𝑡 -1	19.9 m 2 g 𝑃𝑡 -1
	ECSA loss	59 %	59 %

in the Supplementary Material). As shown in Figure

3

, these two techniques indicate an important increase of the nanoparticle size and coherent domain size, respectively. Rietveld refinement of the WAXS patterns gives a coherent domain size increasing from 3.7 nm at BoT to 9.7 nm at EoT (Figure

3

.a), which can be likely explained by both particles aggregation and contact areas with the carbon support (included in the WAXS values but not necessarily electrochemically accessible).

ACKNOWLEDGEMENTS

We would like to thank ESRF for the provision of the beamtime at ID31 beamline as well as Helena Isern and Florian Russello for the help with the beamline preparation. R.C gratefully acknowledge financial support from the French National Research Agency through the HOLYCAT project (grant number n• ANR-22-CE05-0007) and Dr. Bernard Fraisse from the analysis and characterisation platform of the Pôle Chimie Balard. L.C., , K.K, V. M, F.M. and L.D. gratefully acknowledge financial support from the French National Research Agency through the BRIDGE project (grant number n•ANR-19-ENER-0008-01).

AUTHOR INFORMATION