1 Introduction
Located in Central Mexico and in the vicinity of highly populated urban centers, Popocatépetl is an active stratovolcano only 50 km SE from the Mexico City Metropolitan Area (MCMA) and 30 km SW from Puebla. Large pyroclastic flows have reached the MCMA in former epochs (Siebe et al., 1996) and after 70 years of dormancy the volcano reawakened in 1994, starting its current eruption. The current activity is characterized by permanent passive degassing, interrupted by periodic dacitic dome growth episodes accompanied by strombolian-type activity and Vulcanian explosions.
Popocatépetl is known to have one of the largest contributions of volcanic gas emissions worldwide (Gerlach, 1991), with a total SO2-emission between 1993 and 2001 that is estimated to surpass the accumulated SO2-emissions from the climatic 1991 eruption of Pinatubo (Delgado et al., 2001). In current years, Popocatépetl has been estimated to be among the top ten emitters, with around 2.8% of the total global volcanic SO2 emissions (McLinden et al., 2016). The SO2 plume of Popocatépetl has occasionally been detected by air quality monitoring stations in Mexico City (Raga et al., 1999; de Foy et al., 2009) and Puebla (Juarez et al., 2005) and, under certain meteorological conditions, it can impact regional air quality. Ash emissions from the volcano also affect local aviation.
Ground-based ultraviolet (UV) remote sensing techniques have shown to be effective in the measurement of SO2 gas emissions from this active volcano from a safe distance Fickel and Delgado Granados, (2017); Platt et al. (2018); Galle et al. (2010); Delgado et al. (2001); Grutter et al. (2008); Campion et al. (2012); Campion et al. (2018); Schiavo et al. (2019). Additionally, other volcanic trace gases such as HCl, HF, SiF4 have been detected in the Popocatépetl plume by passive infrared (IR) spectroscopy and give valuable information about this volcano (Love et al., 1998; 2000; Goff et al., 2001; Stremme et al., 2011; 2012; Taquet et al., 2017; 2019).
One of the most abundant volcanic gases, carbon dioxide (CO2) is difficult to measure with remote sensing techniques because the high atmospheric CO2 background concentration means that the relative contribution of volcanic plumes to the measured total atmospheric column is generally very small. On the other hand, in situ measurements are often associated with risk and are generally temporally sparse and spatially localized. Nevertheless, the CO2 signal from the Popocatépetl plume was detected in the past by passive IR spectroscopy (Goff et al., 2001) during extraordinarily strong CO2 emission episodes of 1998, for which the high CO2/SO2 ratios were interpreted to reflect a possible episodic assimilation of limestone during magma migration (Goff et al., 2001). These events were associated with very variable CO2/SO2 ratios and therefore the reported CO2 emission rate is most likely not a representative mean value. Emission rates of up to 100 Gg/d (36.5 Tg/yr) were calculated during this exhalation (Goff et al., 2001).
Aiuppa et al. (2019) combined the estimated CO2/SO2 molecular ratio of 8.2 (Aiuppa et al., 2017), which was based on earlier successful measurements Goff et al. (2001), together with the satellite-based SO2 flux estimate by McLinden et al. (2016) and Carn et al. (2017) in order to estimate an average CO2 emission from Popocatépetl of 9.284 Gg/day, 3.4 Tg/yr or 107 kg/s. However, the assumed volcanic CO2 emissions from earlier measurements of the CO2/SO2 ratio and the extrapolation to quiescent average out-gassing phases using more recent SO2 emission estimates could overestimate the CO2 contribution of Popocatépetl to the global carbon cycle Werner et al. (2019).
It is important to measurements (Goff et al., 2001), together CO2 during various days in the typical and more frequent passive degassing state of the volcanic activity in order to obtain more statistically solid results and gain knowledge of the role these prolonged volcanic emissions play in the local, regional and global carbon budget. In recent years, the remote sensing technique using solar absorption high-resolution spectroscopy in the near-infrared has gained sufficient precision to detect enhancements of the CO2 total column of less than one percent (Wunch et al., 2011). Instruments with lower spectral resolution (Gisi et al., 2012) using a similar measurement configuration have been successful in detecting CO2 enhancements in the volcanic plume of Mount Etna (Butz et al., 2017) during field campaigns of short duration. In that study, the CO2/HCl ratios and those of other molecular ratios showed high temporal variability.
In this work, we measured the volcanic gas composition from a fixed site during an extended time period using a high-resolution Fourier Transform Infrared (FTIR) spectrometer that has been used to document not only the atmospheric variability of CO2 (Baylon et al., 2017), but also ozone (Plaza-Medina et al., 2017) and various gases of volcanic origin (Taquet et al., 2019). This station, the Altzomoni Atmospheric Observatory, forms part of the Network for Detection of Atmospheric Composition Change (NDACC) (De Mazière et al., 2018). In this paper we include information about the site, instrumentation (Section 1.1), and about our measurement strategy (Section 2). Further on, we describe how the volcanic gas emission rate is reconstructed from the given information (Section 3.1). For a particular day (26 April 2015), we present the retrieval strategy (Section 2.2) and how to deal with the airmass dependence (Section 2.2), how to improve the precision of the derived molecular ratio (Section 2.3) and calculate the cross-section (Section 3.2) and emission rate (Section 3.3). In that section, we also describe how we combine the CO2/HCl molecular ratios and the HCl emission rates to obtain a statistically more significant estimation of the total CO2 emission (Section 3.5), to provide then a discussion of the principal uncertainties (Section 3.4). The results are finally presented in Section 3.6, discussed in Section 4, and the conclusions are provided in Section 5.
1.1 Site, instrument, and measurements
Our measurements are taken from the Altzomoni Atmospheric Observatory (19.1187°N, 98.655°W, 3,985 m a.s.l.) using a high-resolution FTIR spectrometer from Bruker, model IFS HR120/5, that contributes to NDACC (De Mazière et al., 2018). The high-altitude station is located 11 km away from the crater of Popocatépetl, as can be seen in Figure 1. A solar tracker (Gisi et al., 2011) directs the radiation of the Sun to the interferometer, which measures spectra with a 0.02 cm−1 resolution (optical path difference = 45 cm), as commonly used in the Total Carbon Column Observing Network (TCCON) (Wunch et al., 2011). In the near-infrared spectral region (NIR), an aperture of 0.8 mm is chosen, which results together with the 418 mm focal length of the focusing and collimating mirror in an external and internal field of view of 2 mrad (approx 0.1°). Thus, the field of view diameter is less than 1/5 that of the Sun (the Sun has a size of approx 10 mrad or 0.5°.). The internal divergence and resulting self-apodization are not critical for the used spectral resolution. The remotely-operated FTIR instrument is started manually by an operator and then programmed to continuously measure a sequence of filters and detectors covering the range from 600 cm−1 (17 μm) to 10,000 cm−1 (1 μm). For these measurements in the NIR region, both CaF2 and KBr beam splitters are utilized and an InGaAs detector records the signal of the interferograms (4,000 - 10,000 cm−1). Unlike the MCT and InSb detectors used for the mid-infrared (MIR), the InGaAs detector is not cooled by liquid nitrogen (LN2). All applied detectors are able to record the interferograms in direct current mode making it possible to monitor the solar intensity variations and correct them before the Fourier transformation is performed (Keppel-Aleks et al., 2007). These intensity variations originate from thin clouds, or specifically in our case, from the volcanic plumes that move across the field of view during a single measurement (time of the forward and backward scan is 30-40 s). The correction, which removes smoothed intensity fluctuations of the interferogram outside of the center burst, located at zero path difference, is realized in this work using routines from the “CALPY” software package (Kiel et al., 2016) and slightly improve the precision, Supplementary Section S2.
Measurements are possible almost every morning, but rising clouds typically prevent measurements in the afternoon. As shown in Figure 2, predominant wind directions at pressure level 500 hPa (around 6 km a.s.l.) are west to east and east to west. However for the detection of the plume, wind towards the north (WD=180°) is necessary. Figure 2 shows the wind speed and direction on all days with concurrent HCl measurements. The wind data (500 hPa) are taken for the time at which the highest by PROFFIT96 calculated HCl vertical column of the corresponding day was measured. Of all measurement days, 14% (approx once per week) show a maximum HCl column greater than 1.5E17 Molec/cm2. This value is twice the median value μ (7.6E16 Molec/cm2) and ad hoc chosen to classify days with a clear volcanic signal in the measurements. The corresponding wind direction frequency distribution for days on which a volcanic plume was detected shows a distribution around the wind directions towards North. Applying different, stricter criteria for the classifying of days with volcanic event detection, as 3 or 4 times of the median value results in a subset of 10% or only 6% of all days being chosen, but the relative distribution of wind directions shows a similar pattern and is valid for the subset of 25 measurement days with a HCl slant columns > 1E18 Molec/cm2 which have been chosen for the analysis of the CO2/HCl ratio Section 3.5. The PROFFIT96 output is the vertical column assuming a horizontal homogeneous atmosphere. Therefore it is simpler to classify the days using the vertical columns. However, the signal in the absorption spectrum is proportional to the slant column, which indicates whether we could expect to find volcanic CO2 in the spectra of this day. Most plume intersections are actually recorded by high solar zenith angles (sza > 60°) and therefore the slant column is mostly greater than twice the vertical column, please see Supplementary Material.
2 Methodology
We present a new methodology to determine CO2 emission rates from volcanoes using solar absorption IR spectral measurements taken from a fixed site. Solar IR spectra are continuously measured while the direct sunlight passes through the volcanic plume, allowing for an enhancement in the slant column density of CO2 to be measured as a function of time. Accurate raytracing allowed for retrieving slant columns even for high (>80°) solar zenith angles during sunrise due to the location of the measurement site and the direction (Gisi et al., 2011). Simultaneous detection of other volcanic tracers such as HCl and HF from the same spectra proves that the measured ΔCO2 anomaly is of volcanic origin.
2.1 Solar occultation measurements of gas plumes: vertical column, slant column, plume cross-section and emission rates
In solar-tracking FTIR measurements all detected photons take the same light path and the slant column is well defined by the integration of the gas concentration along the line of sight. The vertical column is the vertical integral of the gas density independent from the observer geometry. The airmass factor is defined by the quotient of slant and vertical column and therefore approximately 1cos(sza) if the curvature of the light path can be neglected, and the atmospheric trace gas distribution is horizontally homogeneous (as is approximately true for background CO2). Vertical columns describe the atmosphere independently of the observer viewing direction and are therefore useful for comparison, validation, bias corrections, or calculation of anomalies. The retrieval code PROFFIT9.6 therefore reports all retrieval results as vertical columns. If a volcanic plume is measured, however, horizontal homogeneity is not a valid assumption and the slant column itself is useful. We can formulate a potential enhancement of the slant column Δcol (sza) by gases in a volcanic plume in the framework of the output of the vertical columns by PROFFIT96.
According to the equation of continuity and the integrated form of the law of Gauss, the emission from a point source can be calculated by the integration of the flux through a closed surface around the point source. Typically we think in a cylinder or cone where only the sides contribute as the top is outside of the atmosphere.
Where E represents the emission rate in molecules per time, A represents the closed surface, ρ is the number density of the considered gas, j⃗ is the flux through the surface, d⃗a is the indefinite small area of the surface. The integration over the area is then divided by the integration along the viewing direction (dz⃗) of an instrument and perpendicular to the viewing direction ds⃗.
Platt and Stutz (2008); Galle et al. (2010); Platt et al. (2018) discuss in more detail different measurement geometries and measurement spectroscopic techniques using the photons scattered or emitted by the sky to calculate these integrals, plume cross-sections, and emission rates, as, for example, i) zenith sky DOAS from mobile platforms or ii) multi axis DOAS (MAXDOAS) using different viewing angles as done by the NOVAC network Galle et al. (2010). For solar absorption spectroscopy in the IR, our viewing angle has to follow the Sun and is identical to the solar zenith angle, and the plume cross-section is given by the integral iii) in Equation 4. It is calculated from the fixed site looking to the Sun and only under favorable conditions once a day. In contrast to the MAXDOAS, we are not able to choose a scanning geometry, but in principle, the plume cross-section is measured in principle the same way. Only a small area contributes to the integral of the closed surface so just this small part of the closed surface has to be measured. For gases with non-negligible atmospheric background concentration such as CO2, the part of the slant column which represents the recently emitted gas is calculated by subtracting the background column. The vertical columns are subtracted and afterward the corresponding slant column anomaly of the target gas is calculated. Inside the retrieval code (PROFFIT9.6), the ray-tracing in the atmosphere is calculated with high precision as already mentioned, but the calculation of the slant column anomaly from vertical column anomaly is done in this work using the simplified airmass factor 1/cos (sza). The solar zenith angle (sza) is the apparent solar zenith angle, the angle at which the instrument is aimed to align the field of view to the center of the Sun and the angle of the line of sight in the lower atmosphere.
The slant columns (HCl) or slant column enhancements above background (CO2) Δcol represent the integrated gas concentration (or gas concentration anomaly) along the line of sight. v⊥ is the wind velocity component perpendicular to the viewing angle and to d⃗s, which corresponds to the connections of two points, each one represents the center of mass of the volcanic gas concentration along the two different lines of sight of consecutive measurements. As we measure the slant column as a function of the viewing and solar zenith angle col (sza) we still have to estimate the distance from the instrument to the plume center to calculate ds (ds = dist ⋅ dsza) between two consecutive measurements and, more specifically the components perpendicular to the line of sight and wind velocity. The wind velocity v⃗ might not be constant with altitude and along the line of sight, but we approximate the wind velocity by querying the NARR model at the altitude of the measured plume center.
Similarly, we estimate the distance between the instrument and the plume center dist. During several minutes required to measure the full plume cross-section, the distance changes slightly as a function of the sza and the spatial extent of the plume. The estimation of distance and height are not independent and a consistency check is possible. The spatial extent and geometrical form of the plume influence the estimation of ds (sza) and therefore the plume cross-section. However, if a constant distance is assumed to be valid for the whole plume cross-section, as is done here, knowledge about the geometrical form of the plume is not needed (Platt et al., 2018).
2.2 Spectral analysis and retrieval strategy
The retrieval strategies used in this work are adapted from Baylon et al. (2017), Taquet et al. (2019) and the standard NDACC and TCCON procedures, thus in this section we only detail the changes from the standard retrievals.
2.2.1 CO2
The analysis of the spectra was performed with the retrieval code PROFFIT9.6 using the CO2 band around 1.6 μm which is used by the Total Carbon Column Observing Network (TCCON) (Wunch et al., 2011), reaching precisions of 0.1% or better in the total CO2 columns (Baylon et al., 2017).
There are some adjustments to the retrieval settings for the volcanic application to the common PROFFIT retrievals at the Altzomoni site. The a priori profile has a constant volume mixing ratio (VMR) of 400 ppm and the chosen constraint is a block Tikhonov type (Rodgers, 2000; Clarmann and Grabowski, 2007), similarly used for retrievals of carbon monoxide within the boundary layer of Mexico City (Stremme et al., 2009; Stremme et al., 2011) but with a free block for the altitudes 4-8 km. This approach ensures that the total column sensitivity, which is also called averaging kernel (AK) is around 1.0 in the lower troposphere and in the altitude of the volcanic plume. The strategy is similar to previous studies conducting either a plume model (Taquet et al., 2017) or scaling the profile of the troposphere (Butz et al., 2017).
2.2.2 O2
In the TCCON and COCCON networks, the vertical O2 column is retrieved from the same spectra as CO2 and used to ensure the quality of the spectrum and to calculate the dry air mol fraction XCO2. In Section 2.3 we describe how we use the O2 column in a similar way to improve the estimation of the CO2/HCl ratio. We slightly adapted the O2 retrieval described by Baylon et al. (2017), in which tropospheric HF is not simulated in the spectra. However, HF is an important interference species and tropospheric HF should be fitted for volcanic application. A constructed HF profile which is constant from 4-7 km (possible altitudes for the volcanic plume) and zero above is scaled in the O2 retrieval. Consequently the interference error between volcanic and tropospheric HF (first-overtone) and O2 is small and irrelevant and no artificial correlation between O2 with any volcanic trace gas (HCl, HF) can be observed.
2.2.3 HCl
HCl is the fourth most abundant volcanic gas emitted by Popocatépetl after H2O, CO2, and SO2 (Love et al., 1998; Goff et al., 2001). Because HCl has a higher solubility in silicate melt than CO2 and SO2, HCl outgassing generally occurs at shallower depths/lower pressure (Aiuppa et al., 2009). In the mid-IR it is normally retrieved at around 2,800 cm−1 together with SO2 at 2,500 cm−1 (Francis et al., 1998; Love et al., 1998; Stremme et al., 2011; Taquet et al., 2019), but can also be retrieved from its first overtone at 5,700 cm−1 (Butz et al., 2017). For this study, we retrieved HCl using the following three microwindows; 5,738.30-5,740.00 cm−1, 5,767.00-5,767.75 cm−1 and 5,779.25-5,779.80 cm−1 with H2O and CH4 as interference gases and the HITRAN 2004, 2006, and 2012 spectral line lists, respectively. The line positions of HCl, H2O, and CH4 have been experimentally slightly adjusted to improve the fit. The HCl retrieval from the first-overtone results in slightly lower column amounts than when using the fundamental band. However, we restricted our analysis to the use of HCl results from retrievals in the NIR region in this work. As a priori, an average of HCl profiles, modeled with the Whole Atmospheric Community Climate Model (WACCM) is calculated and the constraint fixes the retrieved profile above 8 km a.s.l. to this apriori, but allows for free scaling of the lower part from 4 to 8 km. As a consequence, the averaging kernel for total column HCl is close to 1 around the plume height (6 km).
The most important PROFFIT retrieval settings of CO2, O2, and HCl are summarized in Table 1 and more details about the HCl retrieval strategy can be found in Taquet et al. (2019).
2.2.4 Airmass dependence of CO2
It is known that a slightly incorrect simulation of the CO2 line shape and errors in the spectroscopic parameters result in a small artificial airmass dependence of the retrieved CO2 columns. This artifact is especially important and visible at high solar zenith angles (sza) (Wunch et al., 2010). If the chosen a priori CO2 profile differs from the true background atmosphere, the airmass dependence of the sensitivity of CO2 impacts the retrieval as well. Thus the dependence on the solar zenith angle can be eliminated by the use of a specially reconstructed a priori profile ρo (see Supplementary Material). Such correction was applied for the event of 26 April 2015, which is the most relevant CO2 anomaly of our time series, so that it minimizes the air mass dependence and allows for better quantification of the anomaly. As shown in Figure 3, no significant sza-dependance was observed after applying the correction.
To analyze the CO2/HCl ratios in all identified volcanic events, we did not reconstruct an a priori for each individual event. Instead, we fit and subtracted a polynomial of order 2, jointly with the time series of the volcanic tracer HCl and the O2—anomaly and explain the CO2 time series as the instrument viewing direction passes through the volcanic plume as a linear combination of different contributions (Section 3.5). For each measurement, the random and systematic errors in the CO2 column are estimated by the internal error analysis of PROFFIT9.6 (i.e., Barthlott et al. (2015); Plaza-Medina et al. (2017); Baylon et al. (2017)).
2.3 Three strategies to quantify the volcanic CO2 emissions and CO2/HCl ratio
Three slightly different strategies for reconstruction of the CO2 emission rate were applied to data recorded on 26 of April 2015, the day on which the strongest volcanic signal was measured. Strategy i), calculates the plume cross-section directly from the slant column anomalies of CO2 and does not need the retrieval of another volcanic gas. This approach works for the quantification of the CO2 emission even if the composition of the volcanic gas is variable as recently observed by Butz et al. (2017) for Etna volcano.
Strategy ii) assumes a constant molecular ratio between HCl and CO2 in the volcanic gas plume. It calculates the emission flux for HCl, for which the atmospheric background can be neglected, and then multiplies the emission flux by the CO2/HCl ratio. The ratio is obtained as the slope of the fitted straight line in a CO2 HCl correlation plot.
Strategy i) and ii) work only if the CO2 atmospheric background is constant in time. This assumption is rather likely for the remote site Altzomoni and the duration of only a few minutes in which the spectra with a volcanic plume are recorded. Unfortunately, there is a solar zenith angle-dependent bias observed, which is also called air mass dependence (Wunch et al., 2011) as already mentioned above. This solar zenith angle-dependent bias in the CO2 measurements, if not corrected, would affect both strategies (i and ii) as it complicates the determination of the constant atmospheric background of CO2 in strategy (i). This is the reason why an HCl-CO2 scatter plot does not show a very nice and simple linear correlation, which is however needed for strategy (ii).
Knowing the total vertical column of CO2 before and after the field of view intersects the plume allows us to determine the airmass dependence and correct it in the first step. Afterward, in a second step, either the plume cross-section can be calculated from the continuously observed slant columns (strategy i), or the assumed constant CO2/HCl ratio is determined from the CO2 vs HCl scatter plot (strategy ii).
However, the analysis in two steps required by these two strategies can increase the CO2 emission uncertainties and the error propagation from the bias correction to the CO2/HCl ratio is complicated.
Therefore strategy (iii), was designed to retrieve the CO2/HCl ratios from all 25 plume intersections with slant column of HCl > 1018 molecules/cm2. Therefore, the time series of the retrieved vertical CO2 column of each plume intersect is modeled by a linear combination of different contributions from the atmosphere, volcanic plume as well as instrumental and spectroscopic limitations.
1) A atmospheric background CO2 column
2) Solar zenith angle dependence
3) the contribution of the volcanic gas to the measured absorption
4) Differences between calculated and real atmospheric path length
The first two points are modeled by a polynomial. The background column might change either because of a growing mixing layer of Mexico City, which might reach the altitude of the site (Baumgardner et al., 2009; Garcia-Franco et al., 2018), or a different mass of background air might be observed due to the change in the solar angle. In addition, the apparent change in the background column rising from the airmass-dependent bias impacts the behavior of the retrieved background column. A polynomial of the second degree is used to reproduce the function of the solar zenith angle and works for all 25 analyzed events. The pointing error, which affects all gases with an atmospheric background in the same way, can be corrected with the retrieved O2 column anomaly. The TCCON (Wunch et al., 2011) and COCCON (Hase et al., 2015; Alberti et al., 2022) analysis and in the recent work of (Butz et al., 2017) XCO2, the quotient of both columns (CO2 and O2), is calculated and therefore the pointing error is reduced. The retrieval of CO2, O2, and XCO2 time series at Altzomoni and their independent and correlative errors have been studied and described by Baylon et al. (2017). Adding a new parameter (in this case O2 columns) also adds a new error source (Baylon et al., 2017), so that unfortunately a contribution in the error in the XCO2-product arises from the noise in the independent O2 retrieval. Because the cam tracker (Gisi et al., 2011) used at the Altzomoni site is quite accurate, the pointing error is sometimes almost negligible and the use of the retrieved O2 column might sometimes increase the statistical error. An erroneous estimation of the optical path length of the sunlight through the atmosphere affects the retrieved vertical columns of all gases in the same way. Any error in the line of sight results in the anomaly δO2 in the retrieved vertical O2 column, which is expected to be constant in time. This anomaly related to the optical path can be converted to an anomaly δCO2 in the vertical column of CO2 using the factor 400 ppm/21% based on the background concentration of both gases, Eq. 5. We can then subtract the estimated δCO2 vertical column anomaly calculated from the δO2, either completely or also partly (weighting with the parameter α (0 < α < 1). There are plenty of reasons why the pointing error and the measurement noise error vary on some of the 25 events, for example, The impact of the pointing error increases with the solar zenith angle in general, and the signal-to-noise ratio varies depending on the instrument’s optical configuration (KBr or CaF2-beamsplitter), the presence of thin clouds or the volcanic plume might affect the signal to noise ratio. The solar tracking might be degraded through dust and dirt on the optics of the solar tracker and even some transient effects such as oscillations of the solar tracker around the Sun have been observed occasionally. The relative precision in the XCO2 retrieval product might be worse than the CO2 vertical column if the tracker works perfectly, but if the solar tracking is one relevant error source and O2 could be retrieved from spectra with a high signal-to-noise ratio, it is possible to correct for the error in the line of sight successfully. Each short time series of a volcanic event has, according to the different error contributions, one optimal weighting factor α, which will optimize the statistical error in the anomaly of CO2 during the short time that the field of view passes through the volcanic plume.
In Eq. 6 σCO2 is the error in the finalized, optimized and corrected CO2 column, σCO2−independent refers to the statistical error in the retrieved vertical CO2 column, which does not correlate with the retrieved vertical O2 column. σO2−independent is the error in the retrieved O2 column which does not correlate with the CO2 column.
Therefore, we adjust the factor α to minimize the error in the CO2 column anomalies for both higher sza or lower signal-to-noise ratio, which is theoretically given as a function of the correlated error between CO2 and O2 and the independent error of the O2 column. The mean values of the correlative and independent errors have be estimated either by PROFFIT9.6 or the analysis of consecutive measurements in large ensembles (Baylon et al., 2017), but are not known specifically for the measurements of different days and intervals with the volcanic plume measurements.
Figure 3 presents the CO2, O2, and HCl vertical columns, with corresponding error bars calculated with PROFFIT9.6, for the volcanic plume intersection detected on 26 April 2015. The CO2/HCl-molecular ratio was retrieved simultaneously by adjusting of the O2-column anomaly and the sza-dependence as shown in Figure 3. The error in the CO2/HCl-molecular ratio is estimated from the resulting residual as usually done for the measurement noise error in least-squares fitting.
The linear combination of different contributions should approximate the observed time series. The set of factors are the coordinates of the vector x and we describe this linear combination as matrix multiplication of the matrix K and the vector x. The calculated vector x minimizes the cost function CO2(t)−K⋅x2.
The weighting functions (the columns of the matrix K correspond to the background, linear and squared airmass dependence, the proxy for volcanic plume and the proxy for a pointing error in the last column in Eq. 8.
The fourth coordinate from above is the CO2/HCl ratio. With the definition of the Gain-matrix, we can easily calculate the measurement noise error in this ratio as shown in (Eq. 12).
The error in the ratio (Eq. 12) is calculated from the linear regression and helps to weight the events with different strong volcanic signals and quality of measurements.
The CO2 background column can be calculated using the following equation.
Subtracting this background column from the CO2 time series measured on 26 April 2015, we isolate the volcanic CO2 as shown in Figure 4.
3 Results
3.1 The event on 26 April 2015
On the early morning of 26 April 2015, just after sunrise, Popocatépetl was passively degassing as it most frequently does and the sky was clear of clouds. The wind conditions at altitudes between 5,500 and 6,800 m a.s.l. were fairly constant as can be observed from the images captured every minute with the camera at the station pointing towards the volcano (see Supplementary Material). The plume was propagating towards the NE (wind direction = 222°) at a speed of approximately 8.2 m/s according to the radiosonde, launched daily at 06 LT from Mexico City. Continuous measurements of solar absorption spectra resulted in a full cross-section of the volcanic plume at around 06:44 LT, that together with the wind parameters, allowed for the emission flux of volcanic CO2 and HCl to be determined. A schematic of the measurement geometry for this particular day is presented in Figure 1.
3.2 Plume cross-section
As already stated, on the early morning of 26 April 2015, consecutive solar absorption spectra were recorded just after sunrise starting at 06:30 LT. Each measurement had a duration of 37 s. At sunrise, the change in sza is rather fast and a complete cross-section through the volcanic plume was realized in about 15 min. The cross-section is directly given as a set of slant columns for which the time and the sza are known (see Figure 4). The plume cross-section can be converted to a function of height, as shown in Figure 4), since the distance from the observation site to the volcanic plume could be determined based on the known solar azimuth angle, the average wind direction, and the known distance between Popocatépetl’s crater and the observation site (Figure 1).
With the above constraints, we have calculated a distance to the plume of 15.88 km, assuming that the height of the plume is 2.76 km above the observation site and 1.26 km above the crater. That the plume which propagated 18.11 km during about half an hour to reach the instrument line-of-sight, was rising to a height corresponding to its potential temperature (altitude of 6.75 km a.s.l). This propagation can be confirmed by a set of images taken from the Altzomoni Site (See animation in Supplementary Material).
Alternatively, data from the North American Regional Re-analysis Data (NARR) (Mesinger et al., 2006) were used for the determination of wind velocity, plume height, and plume distance, leading to slightly different results (Table 2). This difference is used later to estimate the error associated with wind speed and direction uncertainties. However, the plume was observed at different times and the time intervals between our observation and the last available daily radiosonde observation could be quite large on some days. Instead, we decided to use the NARR data, which are available every 3 h and also assimilate the information of the radiosonde measurements. The statistical analysis that we present is therefore based on the combination of the FTIR cross-section measurements and wind information from model reanalysis data.
3.3 CO2 emission rate calculation
For the calculation of the flux, or the instantaneous emission rate for a particular event, the wind velocity perpendicular to the plume cross-section is needed. This perpendicular component of the wind vector is estimated from the geometry as shown in Figure 1 and based on the radiosonde and/or the NARR model data. For the 26 April 2015 case study, the wind speed and direction used for the calculation are 8.2 m/s and 222°, respectively, leading to an angle between the measurement line-of-sight and the direction in which the plume propagates of 37.4°. The velocity component perpendicular to the cross-section is 8.2/sin (37.4°) = 4.98 m/s. The amount of molecules in the plume slice measured as the Sun increased its elevation angle across the plume is obtained by integrating the SCD over height. The instantaneous emission rate is then calculated by multiplying these two values, thus ECO2 = 36.7 × 1025 Molec/m × 4.98 m/s = 3.0 kMol/s, or in terms of mass, ECO2 = 134 kg/s (116 kg/s using NARR).
3.4 Errors and uncertainties
3.4.1 Gas ratios versus direct CO2 emissions
The main uncertainty is given by the precision in the total CO2 slant column measurement. The volcanic anomaly in the slant column of CO2 is less than 1% of the atmospheric background (vertical column) and the precision of the measurement determined as the random error is around 0.2% of the total slant column, which results in a relative error greater than 20% in each slant column measurement. Figure 4 shows the retrieved slant column anomaly ΔCO2(t)=1cos(sza(t))⋅VCO2(t)−VCO2bkg(t). The errors in this slant column anomaly after background correction and subtraction are estimated as the propagation of the random error (errCO2ST) of the vertical column of CO2 into the slant column anomaly:
The statistical error in the retrieval errCO2ST is analytically calculated by PROFFIT9.6, and the estimated uncertainty in the input parameters and the contributions of different error sources are described in more detail by Baylon et al. (2017). For the full CO2 cross-section the weighted Gaussian sum of all errors in the individual column measurements results in an error of 35%. This large uncertainty can be minimized if we assume that the CO2/HCl ratio is constant during the measurement of a single cross-section. Thus, we can independently calculate the HCl cross-section, which has a better signal-to-noise ratio, and determine a value of CO2/HCl for that measurement sequence. As the atmospheric background of HCl is lower than the detection limit, the HCl slant column amount and the resulting cross-section can be determined with a far smaller uncertainty, and the main contributing error source will thus arise from the molecular ratio. In the slope fitted for the case study, shown in Figure 4 (inset), a relative error of 13% = 1.25/9.61 can be calculated from the 95% confidence interval. The same analysis could be done with HF instead of HCl.
3.4.2 Wind data
The uncertainty in the wind direction affects the emission rate in two different ways: a) The wind direction is an essential ingredient for reconstructing the measurement geometry, as in Figure 1, and impacts the calculation of the distance between the observation site and the plume cross-section. b) The wind direction is used to calculate the component of the wind velocity perpendicular to the plume cross-section as it is needed for the application of Gauss’s Law. In a) the average wind direction of the plume is required, while in b) the wind direction just during the measurement of the plume cross-section is needed.
We assume both the mean and the final wind direction to be 222°, but we count them as two independent components for the calculation of the emission rate. The error in both the mean and the final wind direction is estimated to be ±5°. This ad hoc estimated uncertainty in the wind direction corresponds to the median of the deviations of the two wind directions at 500 hPa and at 450 hPa pressure level of the NARR-reanalysis data. The +5° and −5° error propagates to an error of around +2.9 km or −2.3 km, respectively, resulting in an average of 17% for the reconstructed distance between the measurement site and the plume cross-section. The +5° or −5° error in the final wind direction results in an error in the wind component perpendicular to the cross-section of −11.8% and +11.0%.
Wind speed is critical for the flux calculation. In our case study, the wind speed measured by the radiosonde was almost constant 8.2 m/s over the various altitudes around the plume height (Table.2). The wind speed from the NARR model can be obtained from both its meridional (v) and zonal (u) components as WS = u2+v2 and was found to be slightly lower, 6.8 m/s. This discrepancy of 1.4 m/s was used to estimate the error in the wind speed. For the 26 April 2015 event, a more confident result is expected from using the velocity from the radiosonde since it was taken in the same hour. The total error in the emission rate is thus estimated to be around 30%. These values are summarized in Table 2.
3.5 Time series: HCl emission rates and CO2/HCl ratios
In total, we analyzed 25 occurrences during which our instrument field of view intersected the volcanic plume and during which HCl and CO2 could be simultaneously detected. All of these plume intersections have HCl slant columns greater than 1.0E18 Molec/cm2 but volcanic CO2 anomalies close to their detection limit. Applying strategy (iii) (see Section 2.3) and using scaling factor applied to the HCl time series yields the CO2/HCl molecular ratio and its standard error for each event.
Figures 5A, B show the variability of the CO2/HCl ratio for each plume intersect in the time series. The correlation of each molecular ratio with the HCl emission rate is shown in Figure 5D). The time series does not show a significant pattern or trend. The histogram of the derived ratios might be explained by a log-normal distribution convolved by a normal distribution, since all ratios should be positive, but are determined by relatively large measurement errors (Figure 5D). The weighted average CO2/HCl ratio is found to be 11.4 ± 4.4 Molec. ratio and the correlation plot (inset d)) shows no significant correlation between the emission strength and the Molec. ratios (R2=0.03).
The Figure 5C shows the time series of the 54 events for which the mean emission rate of HCl could be calculated. The HCl flux is rather variable ranges from 0.007 kg/s to 12.6 kg/s, with a median of 2.7 kg/s and an average and standard deviation of (3.0 ± 3.5) kg/s. From those, the 25 events which were also chosen for the estimation of the mean CO2/HCl molecular ratio are plotted in red and marked with a plus in the same figure.
The use of O2 anomalies for the determination of the volcanic CO2 is not crucial for the case study (26 April 2015), but useful for the other 25 volcanic events, which are used to determine a mean CO2/HCl molecular ratio of 11.4 ± 4.4 and have rather large errorbars. Figure 3B shows scattering in the vertical O2 column of around 0.2%, which is a similar magnitude as the volcanic anomaly in the CO2 column. Interestingly, the CO2/HCl molecular ratio determined for the above average emission rate on 26 April 2015 (8.46-9.6) is similar to the average ratio obtained from the full time series (11.4 ± 4.4, see Table 3).
As can be noticed from Figures 5A, B, despite the larger emission on 26 April 2015 (Figure 5C), the molecular ratio of the event lies very close to the mean value.
3.6 Total emissions
The CO2 signature was only found strong enough to calculate the CO2 plume cross-section with a reliable error (35%) during a single plume intersect, occurring during passive degassing. The maximal slant column anomaly of CO2 reached a value of 1.6% above the vertical background column, which is clearly larger than the detection limit and precision. In general, the CO2/HCl ratio is variable from day to day (Figure 5) and initially, it was not clear whether this single measurement occurred through an extraordinary enhancement of the CO2/HCl ratio as found by (Goff et al., 2001) or due to an enhancement of the overall gas emission. On the event of 26 April 2015 the volcanic emissions reflected a typical behavior of passive degassing, during a dome destruction period (photos, Supplementary Figure S1).
Overall, we found that the CO2/HCl ratio on 26 April 2015 was similar to the mean ratio of CO2/HCl for the entire 2012-2016 measurement period. This hypothesis is supported by the analysis of 25 of 54 volcanic events, for which we found an average CO2/HCl molecular ratio of 11.42 ± 4.44. Thus the CO2/HCl ratio of 9.84 ± 1.26 measured on 26 April 2015 seems to be a typical ratio.
The CO2 annual emission rate for the 2012 to 2016 period is estimated from the average CO2/HCl molecular ratio (11.42 ± 4.44) and the HCl emission flux (3.0 ± 0.3 kg/s) and found to be (41.2± 16.7) kg/s.
Using the HCl-product from this work and the SO2 column retrieval from spectra recorded with the NDACC-filter 3 and 6 and detectors InSb and MCT, (Taquet et al., 2019), but measured within a 5 min coincidence criteria, an average molecule ratio SO2/HCl of 6.6 was determined for the time period of this work (Supplementary Section S4). This ratio allows us a comparison with SO2 emission data taken from the literature. Recent studies reporting SO2 emissions of Popocatepetl based on different techniques are summarized in Table.3. Among these studies are UV-cam measurements (of 29 November 2013 (Campion et al., 2015) and 27 days between November and December 2016 (Schiavo et al., 2019), a former study using data of 1 month from various techniques during the MILAGRO field campaign in March 2006 (Grutter et al., 2008), results from the NOVAC-project (Arellano et al., 2021) and with the space-based SO2 emission estimation covering the period 2007-2014 (McLinden et al., 2016). These studies cover a similar period as ours and help us to validate our new method to calculate volcanic gas emissions. The average emission rates of HCl, SO2, and CO2 in the period from October 2012 to June 2016 yield to (3.0 ± 0.3) kg/s, (35.7 ± 4.2) kg/s, and (41.2 ± 16.7) kg/s respectively (corresponding to the annual emission rates of (0.094± 0.011) Tg/yr, (1.13± 0.12)Tg/yr and (1.30 ± 0.53) Tg/yr).
4 Discussion
The linear combination of a CO2 background column, sza-dependence, pointing error, and volcanic anomaly successfully, gives us the confidence to apply this correction to the 25 plume intersections (days) where the highest HCl-slant columns have been measured and the slant HCl columns were greater than 1018 Molec/cm2. The indirect estimate of the CO2 emission, based on the CO2/HCl ratio (11.4 ± 4.4 molecular ratio) and the HCl flux, avoids statistical bias since even events with CO2 anomalies below the detection limit are taken into account in the determination of the emission rate.
4.1 Comparison of volcanic gas flux measurements with other methods
For the comparison of the volcanic gas flux from Popocatépetl with other measurements, in general, only SO2 flux measurements are available from the literature, thus a comparison has to estimate the average molecular ratio of SO2/HCl during passive degassing of the Popocatépetl which is around 6.6 ± 0.3 for the period 2012 to 2016 and using the HCl retrieval from the first overtone. The slight difference in the SO2/HCl molecular ratio values reported by Stremme et al. (2011); Taquet et al. (2019) (6.8 and 9.1) are due to the spectroscopic inconsistency and differences in the time period. Since SO2 can be readily measured with conventional volcano monitoring techniques (e.g., Mobile DOAS, Scanning DOAS and SO2 camera (Grutter et al., 2008; Campion et al., 2018; Platt et al., 2018; Schiavo et al., 2019; Arellano et al., 2021), it is common practice to determine the emission rate with the X/SO2 ratio, where X represents CO2, HCl, etc. To apply this method we need to calculate the CO2/SO2 ratio which we do by measuring the HCl/SO2 ratio and scaling this with the previously determined CO2/HCl ratio. The resulting CO2/SO2 molecular ratio is 2.06 ± 0.7. The uncertainty of 33% in the average ratio is dominated by the uncertainty in the CO2/HCl ratio and takes both into account i) the natural variability of the ratios and ii) the measurement errors. While the derived average CO2/SO2 ratio is in the range reported by Roberge et al. (2009), it is smaller than the estimation for the average ratio by (Aiuppa et al., 2019), as shown in Table 3 and on the lower range from the values of Goff et al. (2001).
In order to validate our method with previous studies, such as the study by Grutter et al. (2008) that combines various methods using thermal emission and UV-based techniques, or the studies by (Krueger et al., 2013), (Lübcke et al., 2013) and (Campion et al., 2015) using 2D visualisation methods with either a thermal emission spectrometer or a UV camera on particular days and (McLinden et al., 2016) or (Carn et al., 2017) basing their studies on long term satellite measurements, we have converted our HCl-flux results to SO2 fluxes. HCl and SO2 are both volcanic gases with low tropospheric background, however, SO2 is also emitted by anthropogenic sources (de Foy et al., 2009). HCl, on the other hand, is less abundant in volcanic plumes, can only be measured effectively by IR spectroscopy and has not yet been measured by imaging techniques. The HCl-flux in this work multiplied by an average molecular SO2/HCl ratio of 6.6 and corrected for the different molecular weights of SO2 and HCl correspond to a SO2 flux of (35.7 ± 4.2) kg/s which can be compared to different measurements and estimations reported in the literature and summarized in Table 3. These reported SO2 emissions range from 5 kg/s to 579 kg/s in different periods. Especially if we consider the activity of the Popocatepetl during the period 2012-2014 (McLinden et al., 2016; Carn et al., 2017; Boulesteix et al., 2022), we find values corresponding to the 35.7 kg/s obtained from this work. The methodology presented here benefits from the location of both the Altzomoni Atmospheric Observatory and the volcano, and from the typical wind direction during the early morning driving the plume in the line of sight of the instrument and allowing a large number of solar absorption plume measurements. Our HCl-fluxes time series did not capture plumes from any strombolian phases associated with dome growth episodes or Vulcanian explosive activity accompanying dome destruction. In the long-term record both of these episodes are relatively rare compared to passive degassing and excluding these episodes of Strombolian and Vulcanian-style eruptions will bias the estimation only slightly. In and after the explosions a larger amount of gas is released and therefore a positive anomaly is expected, but prior to these explosions a decrease in gas emissions is typically observed. Somehow gas is retained until the explosions and the negative anomaly in the emission just before this explosion compensates for the positive anomaly due to the explosion (Love et al., 2000; Taquet et al., 2017). Our SO2 flux of (36 ± 4) kg/s agrees well with the 31.2 ± 2.4 kg/s average reported for the 2012-2015 period Supplementary Material of Carn et al. (2017) and Boulesteix et al. (2022) and supports the assumption that quiescent emissions dominate the overall emissions from Popocatépetl.
4.2 Quantification of volcanic CO2/HCl ratio by a linear joint fit
Measuring the CO2 flux and the CO2/HCl ratios is challenging because the volcanic part of the slant column of CO2 is in the subpercentage region. The high CO2 background and the sza-dependence for measurements with low solar elevation angles are much larger than the volcanic contribution to the total column as described before. However, a joint fit of different proxies for different contributions can explain the measuredtime series of CO2. Here the HCl time series is used as proxy for the volcanic plume, the columnar O2 anomalies (α400ppm21%δO2) as proxy for the pointing error, which refers to an error in the optical pathlength in the atmosphere and a polynomial covers the remaining contributions as a linear sza and squared solar zenith angle dependence and the background CO2 column. This strategy of fitting and subtracting the O2 anomaly instead of using the XCO2 product is based on the recent work by Baylon et al. (2017), in which the error budget is separated into correlated and independent errors. Calculating the quotient XCO2 is similar to a subtraction using alpha=1.0, as CO2O2+δO2≈CO2O2−CO2−(400ppm/21%)δO2O2. We cannot explicitly explain why α is around 0.25 on the presented case study of 26 of April 2015, and varies between 0 and 0.6. But there are various aspects that impact the retrieved value of α: i) the ratio between the noise error in O2 and the O2-CO2-correlative error, ii) the systematic parts between the O2-CO2-correlative error, iii) differences in the strength of the O2-CO2-correlative error, as maybe raytracing in the 7,800 cm−1 is slightly different from the 6,300 cm−1 spectral region, thus differently sensitive to pointing errors produced by the presence of a volcanic plume and other effects.
During decompression of magma rising from depth towards the surface, CO2 first partitions to the vapor phase, only later followed by sulfur and chlorine (Edmonds et al., 2022). We, therefore, expect that deep magmatic recharge could lead to an increase in the CO2/HCl ratio in the gases emitted from the top of the volcano. On the other hand, efficient release of chlorine from the magma can only occur as it reaches very shallow depths or even atmospheric pressure (Edmonds et al., 2001). Therefore, the HCl flux will primarily be driven by the availability of fresh magma at the surface.
Our measurements provide some limited information on magma migration in Popocatepetl’s plumbing system. For one, the variability of the HCl flux is significant and exceeds the errors associated with each individual measurement. This variability is likely related to very shallow processes, in particular, the volume of chlorine-rich magma at or very close to the surface at the time of measurement. The variability of the measured CO2/HCl ratio, on the other hand, is mostly smaller than the measurement uncertainty. We therefore cannot interpret individual measurements with regard to their underlying magmatic processes. We can, however, establish a baseline for the CO2/HCl ratio to which future measurements can be compared. Based on the range of values we detected (see Figures 5A, B), any future measurements of CO2/HCl significantly higher than 100 (after accounting for measurement accuracy and precision) should be considered anomalous and could be an indicator of deep magmatic recharge at Popocatépetl. Thus, our measurement techniques may provide valuable information on magma migration and possible future eruptive episodes.
The CO2/HCl (11.42) molecular ratio might be converted to a CO2/SO2 ≈ 1.7 using the SO2/HCl ratio of 6.6. Our average CO2/SO2 ratio is consistent with that measured by Roberge et al. (2009) in a melt inclusion (i.e.,: ≈1 Molec. ratio for an inclusion captured at 400 MPa (15 km) representative of primary magma). The CO2/SO2 molecular ratio of around 2 is however in general lower than those obtained in studies carried out in 1993 and 1994 using remote sensing (Goff et al., 2001). Since we also see a large variability in our results and a distribution with a larger tail towards positive anomalies (Figure 5B), this ratio might have significantly higher values during some periods of short duration. However, in this work the number of measurements of the ratio is limited and the precision of the retrieved ratio does not allow to determine the variability, just the overall background value of the CO2/HCl ratio in the plume of Popocatépetl.
5 Conclusion
In this contribution, we determine the average CO2 and HCl emission rates of the active Popocatépetl stratovolcano by near-infrared solar absorption spectroscopy from a fixed site. Our measurements indicate that Popocatépetl has emitted in recent years an average of 1.3 Tg of CO2 and 95 kt of HCl per year. The CO2/HCl molecular ratio is on average around 11.42 but shows a large variability. No significant correlation between the HCl flux and the CO2/HCl ratio was found. A particular case study is presented for a measurement recorded during the early morning of 26 April 2015. An anomaly in the CO2 column on this day was visualized and the emission rate of Popocatépetl could be calculated directly from the CO2 column measurements. Three different strategies for CO2 and HCl flux analysis have been compared and cross validated. For the estimation of a representative mean emission rate for HCl and CO2 during the period between 2012 and 2016, only one of these strategies were followed since the CO2 signal was not always strong enough for a direct estimation. During 54 days in the measurement period, a cross-section of the volcanic plume was measured and the HCl emission rate was determined. The CO2/HCl ratio could be measured in about half of these events and a mean emission rate of CO2 could be calculated using this ratio and the HCl-flux.
Assuming a SO2/HCl molecular ratio of 6.6, the average HCl emission during 2012-2016 obtained in this study corresponds to an annual emission rate of (36 ± 4) kg/s (3.1 ± 0.4) kt/d or (1.14 ± 0.13) Tg/yr) of SO2. This result is completely independent of other methods used so far and is in agreement with other reported SO2 emission rates. The CO2/SO2 molecular ratio of 1.7 is in accordance with some ratios measured by Roberge et al. (2009). However, the 41.2 kg/s average CO2 emissions obtained in this study (equivalent to 1.3 ± 0.5 Tg/yr) is significantly lower than what has been estimated and assumed in former studies (Delgado et al., 2001; Goff et al., 2001; Witter et al., 2005; Aiuppa et al., 2017). According to these new results, the total CO2 emissions of Popocatépetl contribute only with about 1% and 3.4% of the estimated anthropogenic CO2 emissions in the Mexican territory (130 Tg/yr) and in the Mexico City Metropolitan Area (44 Tg/yr). Nevertheless, these CO2 emissions from volcanic origin can be as strong as a coal-fired power plant and should be considered in regional carbon cycle studies and also included in global volcanic emission inventories.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author contributions
WS coordinated the work and carried out the early morning. MG contributed in the reconstruction of the volcanic emission rate from solar absorption measurements, combining perpendicular cross-section, the geometry and writing of the manuscript. JB retrieved the CO2 and O2 columns from the NIR spectra. TB leads the German-Mexican collaboration for long-term operation as part of NDACC. FH has helped in the design and setup of the spectrometer and solar tracker and has developed the retrieval code PROFFIT. NT is responsible for the volcanological aspects and interpretation of the gas composition, reviewing the retrievals and the results. She has contributed to the analysis and writing parts of the paper. AB is responsible for the data acquisition, managing and processing of the spectra, as well as in the operation of the spectrometer. BS reviewed the emission of Popocatepetl and calculated the SO2 emissions. EP-M helped to set up the spectrometer and carried out most of the solar FTIR measurements in Altzomoni. She contributed to the analysis, interpretation and writing of the paper. CR helped with the estimation of the SO2 flux from Popcatepetl. CR has contributed to the analysis, interpretation and writing. All authors contributed to the article and approved the submitted version.
Funding
We acknowledge the financial support of DGAPA-UNAM PAPIIT (IN107417, IN111521, IA101620) and CONACYT (275239, 290589).
Acknowledgments
WS thanks DGAPA PASPA for supporting this research during the sabbatical 2019/2020 and The IMK-KIT for hospitality. CONACYT is acknowledged for the studentship of JB. We are grateful for the technical assistance provided by Alfredo Rodriguez, Miguel Robles, Delibes Flores, Omar López, and the Insituto de Ciencias de la Atmósfera and Cambio Climatico (UNAM) for the institutional support to carry out this study. Especially we thank Dr. Thomas Boulesteix for his help at the Altzomoni site, his fruitful discussions, and his motivation. We thank both reviewers for their constructive critics and helpful comments. Especially we thank Nicole Bobrowski (reviewer 2) for her good ideas and very careful revisions which have improved the manuscript substantially.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The reviewer AB declared a past co-authorship with the authors FH, TB to the handling editor.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1022976/full#supplementary-material
