ORIGIN AND DISTRIBUTION OF DIFFUSE SOIL CO2 GAS EMISSIONS ACROSS TURRIALBA AND IRAZÚ VOLCANOES, COSTA RICA

............................................................................................................................ vii


Author Contribution Statement
The Soil Carbon Isotope Measurements portion of the Methods section was written by Christofer Jimenez of the Observatorio Vulcanológico y Sismológico de Costa Rica and the National University of Costa Rica. All other sections were completed by the author. vi

Abstract
We characterized large-scale volatile emissions across the summit and flanks of the actively degassing Turrialba volcano, Costa Rica, using soil gas flux measurements and 13 C isotopes. The objectives of this study were the following: 1) to monitor changes in the magmatic activity and identify source contributions over Turrialba using measurements of CO2 soil gas emissions since the volcano last erupted in 2014-2015, and 2) to identify the location and extent of magma at depth, and the structures that allow gas transport to the surface. Degassing at the summit is concentrated along one normal fault lineament, the trace of which runs through the west active crater and fumarolic fields which lay on either side. The high soil gas flux measured, visible evidence of active degassing through fumaroles, and magmatic δ 13 C (up to -1.19‰) suggest permeability along this structure that gas follows from magma at depth. δ 13 C samples collected in and around fumaroles at the summit are more enriched than those collected within an area known for consistent fumarolic activity in the Ariete fault, which has also been used by previous studies as a baseline for hydrothermal-magmatic δ 13 C signatures (~-3‰). Additionally, there is a higher CH4/CO2 ratio present in the Ariete fault than at the summit. The difference in isotopic signatures between flank and summit fumaroles is interpreted to indicate magmatic-dominated activity at the summit of Turrialba, possibly comparable to that present at the end of the 2014-2015 eruptive cycle. A similar zone of high diffuse degassing with a more hydrothermal carbon isotopic signature (~-4‰) occurs in the depression between Turrialba and Irazú volcanoes, presenting as a hot spring known as Santa Teresa. We hypothesize the origin of CO2 flux in this area to be related either to transport through the hydrological system from higher in the volcanic system or from previously unmapped magma at depth below the depression. More detailed mapping of soil CO2 emissions must be done between Turrialba and Irazú in conjunction with geophysical monitoring to further advance our understanding of the magma distribution and potential magmatic connection between these volcanoes.

Introduction
Volcanic systems emit volatiles at varying levels through advective degassing in active plumes or fumaroles and through diffuse soil gas emissions. Carbon dioxide (CO2) is a primary component of this gas, and the carbon isotopic composition has been shown to reflect the origin and character of degassing on active stratovolcanoes. Volcanic systems exsolve CO2 at depth which maintains a specific isotopic signature (relative to nonvolcanic CO2) as it rises and reaches the surface. While the primary surface expression of volcanic gas is as often a volcanic vent plume or fumarole, there is a significant anomalous diffusive component which is often more difficult to constrain, and in some cases even greater than any advective component. There is a well-established link between subsurface structures (i.e. faults, fractures, or collapse structures) and diffuse degassing, known as diffuse degassing structures (DDS) (Cardellini et al., 2017;Chiodini et al., 2001;Chiodini et al., 1998;Lin et al., 2019;Melián et al., 2014;Notsu et al., 2006;Viveiros et al., 2020;Werner et al., 2014;and others). Volcanic CO2 tends to follow DDS rather than being distributed homogeneously across the volcanic flanks, though distributions of diffuse soil gas emissions do not occur at the same rate or evenly along a fault/fracture (Jolie et al., 2016). This behavior has been recorded at Yellowstone volcano, USA (Lin et al., 2019), Campi Flegrei caldera, Italy (Cardellini et al., 2017), and in the East African Rift (Jolie et al., 2019). Epiard et al. (2017) have also shown that diffuse degassing is still a significant contributor to the total carbon budget at Turrialba volcano and others along the Central American Volcanic Arc. Studies have shown that determining changes in the CO2 flux rate and distribution can improve our understanding of concomitant changes in the subvolcanic magma reservoir that might signal volcanic unrest (Camarda et al., 2012;Melián et al., 2014;Seiler et al., 2017). Increased rates of CO2 flux (fCO2), for example, have been correlated with pre-eruptive activity at the El Hierro volcano (Melián et al., 2014) and Etna volcano (Seiler et al., 2017). Additionally, CO2 flux can provide significant constraints on magma location at depth and constraints on when new material may have been injected (Camarda et al., 2012;Notsu et al., 2006;Werner et al., 2014). However, where a biological layer exists between the source of gas emissions and the surface there can be mixing between CO2 of volcanic, biogenic, and atmospheric origin (Chiodini et al., 2008;Malowany et al 2017;Viveiros et al., 2020). Carbon dioxide sources can be roughly distinguished from one another by using a probability plot (Cardellini et al., 2003;Chiodini et al., 1998) or through integrating soil CO2 flux and carbon isotope measurements (Chiodini et al., 2008;Viveiros et al., 2010). Carbon-13 in CO2 (δ 13 C) is well established as a marker of CO2 origin and has been used with great success as an indicator of changes in source conditions of volcanic systems. For example, levels of hydrothermal-magmatic activity (Moussallam et al., 2014), and separation of biogenic influences on the volcanic gas composition (Chiodini et al., 2008;Viveiros et al., 2020). Chiodini et al., (2008) have shown a clear difference between biogenic and hydrothermal CO2 at the Solfatara of Pozzuoli volcano. δ 13 C has also been used to track changes in the volcanic system at Turrialba periodically throughout the last 20 years. Vaselli et al., (2009) record changes in the isotopic signature across the summit from 1998 until 2008, with Moussallam et al., (2014) measuring in the summit vents in 2013 and Malowany et al., (2017) continuing measurements in 2014. These studies have shown δ 13 C to be an effective indicator of the balance between magmatically-dominated activity and a periodic hydrothermal reservoir, when combined with measurements of other gases or elements (i.e. plume SO2, He, and CO2 flux).
Turrialba volcano displays cyclic behavior between hydrothermal-magmatic dominated activity, related to phreatic eruptive activity, and magmatic-dominated activity, related to phreatomagmatic activity. Additionally, De Moor et al. (2016) suggested the existence of two magma reservoirs beneath the volcanic summit, one deep (~8-10 km) and one shallow (~3-5 km). Similarly, at Mt. Etna volcano, Italy, Camarda et al. (2012) showed that peripheral, or flank, contributions to CO2 flux may be indicative of a deeper magmatic source, whereas measurements at the summit or higher on the volcanic edifice may represent a shallower magma reservoir. This scenario may also exist at Turrialba, and if so, measurements of flank degassing may provide insight to the deep magmatic system once origins of CO2 flux are well constrained. At the minimum, diffuse gas flux measurements across Turrialba will provide further constraints on magma chamber geometry and magma migration, realities which may have been over-simplified in the past. However, because Irazú volcano is located relatively nearby, measurements of fCO2 within the depression spanning the two volcanoes may include gas emissions from one of either volcanic reservoir beneath Turrialba volcano, a combination of both, and may include a contribution from the Irazú volcanic reservoir as well.
In this study, we combined measurements of diffuse soil CO2 emissions, soil isotopic carbon, and CO2 and CH4 concentrations, on the Turrialba volcano and through a topographic depression between Turrialba and the Irazú volcano. The objectives of this study were the following: 1) contribute to a characterization of the state of activity within the Turrialba volcanic system, 2) constrain large-scale diffuse degassing over known and unknown DDS both at the summit and on the volcanic flanks, and 3) explore the applications of δ 13 C with CH4 and CO2 as indicators of CO2 gas origin especially in areas where volcanic activity is not easily recognized.

Geologic Setting
Turrialba volcano lies along the Central American Volcanic Arc, created as the Cocos plate is subducted beneath the Caribbean plate at a low angle (Peacock et al., 2005). It is a freestanding andesitic to basaltic-andesitic stratovolcano, located within 10 km of the Irazú volcano, another active stratovolcano of the same arc ( Fig. 1) (GVP 2021). Regional faulting trends NW-SE, parallel to subduction, with an additional fault system running perpendicular. Both Turrialba and Irazú formed as part of a pull-apart structure trending SW-NE, creating space within the crust for magmatism that would eventually form both volcanic structures (Montero, 2011).
There are three craters at the summit of Turrialba volcano: West, Central, and East. Of the three, only the West and Central craters have recorded fumarolic activity (Vaselli et al., 2009). These craters are roughly aligned SW-NE, where present activity is concentrated at the West "active" crater. Between Turrialba and Irazú lies a large depression (~10 km wide) attributed to glacial activity (Calvo et al., 2019). Within this valley lies a series of drainages and cold-and hot-CO2 springs. Each volcano is known to have an individual magmatic and hydrothermal system, however, there is a possibility that Turrialba and Irazú may share some component of their magmatic plumbing system (e.g. Mt. Saint Helens and Mt. Adams) (Hill, 2009) Therefore, this study investigated not only the summit and flanks of Turrialba and Irazú, but the depression spanning the area between these volcanic massifs. Figure 1. a) Turrialba is on the southernmost end of the Central Cordilleran volcanic chain, part of the Central Volcanic Arc, and formed as a result of subduction off the eastern coast of Costa Rica. b) Three craters make up the volcanic summit; the West, Central, and East. Two fumarolic fields lie to the north and south of the West active crater, outlined in grey. Vaselli et al., (2009) collected samples from the central (light purple point), west (dark purple), and west high temperature (yellow point) fumaroles in 1998-2008. c) Turrialba and Irazú volcanoes are separated by a deep depression, visible from the south. The Santa Teresa spring, discussed in detail below, is located within this area. The dark grey box indicates the summit area of Turrialba volcano.

Eruptive History
The most notable eruption of Turrialba volcano in recorded history occurred from 1864-1866, after more than 100 years of quiescence. This eruption was characterized by ashfall that reached the Central valley, accompanied by a few pyroclastic flows and lahars (Gonzalez et al., 2015). Activity in the area resumed relative quiescence until fumarolic activity was observed beginning in 1980, with an increase in degassing activity beginning in 1996 (Martini et al., 2010); indicative of reactivation of the volcanic system likely initially related to injection of a small batch of juvenile magma (Campion et al., 2012). Since then, the dominant system has been cyclic, alternating between hydrothermal and magmatic dominated systems (de Moor et al., 2016;Malowany et al., 2017;Martini et al., 2010;Vaselli et al., 2009). Between 1998 and 2008 gas chemistry has indicated a transition from hydrothermally-dominated activity (1998 to late 2001) to hydrothermal/magmatic activity (late 2001 to 2007) to magmatic dominated (2007-2008) (Vaselli et al., 2009). We infer a recharge of the hydrothermal system between 2008 and early 2010, with steadily increasing volcanic activity culminating in a phreatic eruption on 5 January 2010 (Alvarado, 2021). This eruption is indicative of a magmatic-dominated system, coinciding or creating a new vent on the inner wall of the SW crater and allowing open-system degassing to begin (Campion et al., 2012). A second phreatic eruption occurred at Turrialba in 2012 after liquid sulfur was observed pouring from the Quemada fracture (a 4-m wide fracture between the central and southwestern craters, trending roughly NE) and vent incandescence was observed on the Southern wall of the Southwestern crater (Global Volcanism Program, 2013). The 2010 eruption marks renewed eruptive activity at Turrialba which still persists.
Following a period of quiescence, a period of intense phreatic to phreatomagmatic eruptive activity occurred from October 2014 to May 2015, beginning with a blast strong enough to cause the eastern wall of the west active crater to collapse. Two cycles of degassing are evident over the subsequent eruptive activity and are interpreted to be linked to pulses of new magma toward the summit. A hydrothermal gas overprint is evident in the first eruptive phase, but not in the second, suggesting "boiling off" of the hydrothermal reservoir or a change in magma pathways; marking the transition from phreatic to phreatomagmatic activity (de Moor et al., 2016). Summit vent degassing with periodic ash-fall continues to present as part of the ongoing eruptive cycle at Turrialba.

Soil Diffuse Degassing Measurements
The accumulation chamber method has proven to be an effective tool for measuring diffuse CO2 gas flux over volcanic flanks (Chiodini et al., 1998;Epiard et al., 2017;and others). We followed the methodology well-established by previous studies to measure CO2 flux (fCO2) over 600 points across the flanks of Turrialba and Irazú volcanoes in March of 2021. Two West Systems® flux meters attached to type B accumulation chambers (the larger of two configurations, advantageous for collecting and quantifying high gas flux values) facilitate gas sample flow from the soil-atmosphere boundary within the chamber opening over a LI-COR 820 CO2 gas analyzer (dynamic range from 0 to 20,000 ppm CO2 with an accuracy better than 1.5% of reading). The specifications for each fluxmeter used are comparable. Samples are measured as the concentration of CO2 per second as gas is pumped over a sensor before returning to the sampling chamber (West Systems® Handbook). A "flux", in ppm s -1 , is derived from the slope of a best fit line (R 2 ≥ 0.9) over the course of measurement, which can take >2 minutes at low flux locations. Synchronous measurements of near-surface ground temperature using a Hanna K-type thermocouple at ~15cm depth were collected within 1 m of each point unless ground hardness prevented measurement. Location and elevation values were recorded using a Bad Elf GPS Pro+ where accuracy averaged within 5 m and never exceeded 10 m. Data was collected using the real-time interface between the Bad Elf GPS unit and ESRI ArcGIS Collector offline mapping software. Additional meteorological conditions were recorded, including wind speed and direction using a hand-held Kestrel 3000 wind meter, soil moisture from a Campbell Scientific HydroSense II handheld soil moisture sensor, and barometric pressure and air temperature recorded using internal barometric pressure sensor and thermometer within the Bad Elf GPS, though meteorological conditions remained relatively constant between days of measurement. Control point measurements were taken from the same location at the start and end of each day of field work. These measurements are meant to act as a characterization of variance in flux due to external influence, as we would otherwise expect it to remain relatively stable.
The flanks of Irazú and Turrialba range from agricultural and pastoral fields to natural forest to unvegetated areas of volcanic deposits. Because the study area is so large and heterogeneous relative to a standard study characterizing a single degassing area, a regular geometric measurement grid (which may be preferable for ease of statistical interpolation) is not feasible. Instead, measurement points are guided toward areas of active degassing by known fault maps and the knowledge of local communities, gained through open conversation. In this way, we focus on characterizing areas of high degassing while measuring the spaces between them at a coarse resolution (up to >100 m) as a representation of background CO2 concentrations. Once an area of interest was identified (for example, over high flux areas such as a fumarolic field or actively degassing hot spring), we collected high-density measurements (as fine as 3 m apart) moving outward until flux and concentration values returned to background or near atmospheric levels.

Soil Carbon Isotope Measurements
Samples were collected from the 9 th to 17 th of March 2021 from fumaroles at the summit of Turrialba volcano and from soils on the flanks of Turrialba and Irazú volcanoes. Fumarole gas emissions were sampled by inserting a titanium tube, approximately 2 m in length, into the fumarole to prevent air contamination. Silicon tubing was used to connect the titanium tube to a three-way valve and a 60 ml syringe to extract the gas. The line was cleaned at least three times (extracting and flushing gas with the syringe) before collecting the sample to avoid air contamination. Gas was injected into re-evacuated septum vials (samples of 15ml and 120ml were taken) and overpressurized to ensure successful sample extraction in the lab and minimize air entrainment during transport and storage. At diffuse/soil degassing areas, we collected gas samples by inserting an aluminum tube of approximately 1 m with holes drilled into the lower 30 cm of the tube. Gas was sampled using a 60 ml syringe via a three-way valve. The extracted gas was injected into 15 ml septum vials, after flushing the line three times to avoid air contamination.
Samples were analyzed for carbon isotope compositions of CO2 and CH4 on a Picarro G2201-i isotopic analyzer operated in high precision mode at OVSICORI-UNA (Observatorio Vulcanológico y Sismológico de Costa Rica) following the methods of Malowany et al. (2017). Reliable δ 13 C values can be achieved on this instrument for analytical CO2 concentrations between 400ppm and 2000ppm and for CH4 at 1.5ppm to 20ppm. The samples were extracted from the vials, transferred to 2L flexfoil bags, and diluted using ultra high purity N2. The flexfoil bag was then connected to the Picarro inlet and analyzed for approximately 10 minutes. Standards were run every 8 samples, including NOAA standard tank containing 1.88 ppm of CH4 with ð 13 CCH4 of -47.4 ‰ and 394.8ppm CO2 with δ 13 CO2 = -8.3 ‰ (e.g. Caballo-Chaves et al., 2020) as well as an internal standard tank (100% CO2) with ð 13 CCO2 of -3.8 ‰, which was diluted with UHP N2 for analysis. Concentration and isotopic average and standard deviation values were recorded from the Picarro software using a ~5-minute data window. Sample data were then corrected using the standard data from the same run, and CO2 carbon isotope values are reported relative to Pee Dee Belemnite (PDB). Methane concentrations in the samples were typically too low for reliable determination of ð 13 CCH4, and therefore ð 13 CCH4 values are not reported.

CO2 flux data analysis
Flux values are converted from ppm/s to g/m 2 /day to standardize units. This is accomplished by factoring in barometric pressure, air temperature, and accumulation chamber dimensions following the A.c.K equation outlined in the West Systems® handbook: = 86,400 * 10 6 * * * (1) where P is the barometric pressure in mbar; R is the gas constant of 0.08314510 bar K -1 mol -1 ; Tk is the air temperature in Kelvin; V is the net volume of the chamber, 6.186*10 -3 m 3 for a type B chamber like this one; and A is the chamber inlet area which is 3.140 * 10 -2 m 2 for this type B chamber. Dimensions of the flux in mol m -2 day -1 are found using the following A.c.K.: and from here, units can be converted to g/m 2 /day using the molar mass of carbon dioxide. A value of 1 was added to all resultant flux values for analysis and removed after calculations were completed. Flux values under the detection limit are assigned a value of 1 and remain in the data set as a representation of areas with very little contribution to flux. However, where our data is compared with those of other studies these values are not considered.
Typically, one would interpret the spatial distribution of the diffuse degassing contribution from a statistical analysis of the space between points. However, the spatial continuity is poor due to the inconsistency and irregularity in measurement spacing. The graphical statistical approach (GSA) was used instead; after the methods of Cardellini et al., (2003). The GSA method partitions fCO2 data into groups based on the lognormal probability plot, where each group represents a different composition of gas source biogenic, mixed, or volcanic). This method is advantageous in this case as it does not consider spatial correlation within the data.

Error Analysis
Meteorological variables can potentially influence soil CO2 flux values as they can change over relatively short time frames and impact soil porosity and diffusion between the soil and atmosphere (Camarda et al., 2006;Chiodini et al., 1998;Hinkle et al., 1994). To identify any confounding influence on gas flux measurements, meteorological variables were recorded synchronously with each flux measurement and treated as predicting variables on fCO2 in a multivariate analysis. If flux values were significantly impacted, we would expect a statistically significant correlation to become evident. This analysis is made more robust by employing a machine learning component; especially as there is a high likelihood of adding more data points to the model in the future.

Field Observations
The study area can be divided into six sections based on large-scale topography, land cover type, and volcanic activity (Fig. 2). In areas 1 (Turrialba's summit) and 2 (the Ariete fault or Falla Ariete) gas isotopic composition and CO2 fluxes indicate a strong magmatic source component. The Ariete fault is also expressed as a fault scarp and intermittent stream, with activity concentrated within fumaroles at the topographic low. Area 3 (La Silvia) is defined by a volcanic ash deposit and dead vegetation. This was the site of mass vegetation die-off beginning in 2005, following the initiation of a summit vent plume on Turrialba. Lethal levels of CO2 were present in small, localized areas in addition to volcanic plume gas species (predominantly SO2) carried by persistent trade winds (Global Volcanism Program, 2013;Tortini et al., 2017). Area 4 is characterized by a significant depression where the study area is at least 500 m lower than at any other section and degassing is localized to areas of surface water, including the Santa Teresita hot spring and a series of creeks and rivers. This area is between the Turrialba and Irazú volcanoes. Areas 5 and 6 are defined by predominantly agricultural and pastoral land on the flanks of Turrialba and Irazú. Area 5, known here as "Irazú", is located on the summit and southern flanks of Irazú volcano. We did not measure any significant degassing in this area, allowing any influence from the magmatic system of Irazú to be discounted in this area. Area 6, or "Santa Cruz", is over the Southern and Eastern flanks of Turrialba volcano, excluding the area which is known to contain the actively degassing Ariete fault. Figure 3 displays a visual representation of each area discussed.

Diffuse soil gas fluxes and error
The distribution and range of fCO2 measurement points between study areas is displayed in Table 1. The minimum of any area is below the detection limit of the sensor, and the maximum of any point was 7917.8 g/m 2 /day. Using a lognormal probability plot, we can partition the CO2 flux data into populations based on the most likely origin of the gas (Cardellini et al., 2003). This approach is widely accepted in areas where both magmatic and hydrothermal systems are degassing in conjunction with a biogenic source of CO2. We assume the source of CO2 flux is either magmatic-hydrothermal or biogenic in origin, with the majority of points falling between, where the sources are mixed. Typically, lower fluxes are associated with a biogenic source whereas a higher source is associated with a volcanic origin (Chiodini et al., 1998;Lewicki et al., 2005; and others). Thresholds are determined from a lognormal frequency plot, where flux populations are more easily visualized (Fig. 4). Different populations are defined by bounding inflection points (Sinclair, 1974). Once inflection points are well defined, they can be considered as thresholds. Then, individual measurement can be classified as background, mixed, high, or above the limit of diffuse degassing, therefore classified as advective. Results of this classification are presented in Table 2.  Table 2. Flux point statistics, as determined from a cumulative frequency distribution plot (Fig. 4) after Cardellini et al. (2003) and Sinclair (1974).

Percent of Total Points
Background ≤ 0.5 33.6 % Mixed 0.5 -2.0 60.6 % High 2.0 -4.5 3.5 % Advective 4.5 -9.0 2.3 % 27 measurements of fCO2 and all meteorological variables were collected at a control point in the morning before field work and in the evening afterward (Fig. 5). Values of fCO2 ranged from 0.3 g m -2 day -1 to 5.89 g m -2 day -1 and averaged at a value of 1.89 g m -2 day -1 . Using all fCO2 measurements over Turrialba that included concurrent meteorological measurements of all conditions, a total of 395, as a training dataset we used 12 multivariate regression models to determine if any measured meteorological conditions significantly impacted CO2 flux. The WEKA machine learning graphic user interface was utilized for ease of analysis (Frank et al 2016). Initially, a filter-based attribute selection was used to ensure variable selection was independent of any algorithm. The CFS Subset Evaluation evaluates the worth of a subset of attributes through consideration of predictive ability each secondary variable has on the variable of interest, where more is better, and the degree of redundancy the secondary variables have between each other, where more is worse. The "best" secondary variable is that which has the highest correlation with the variable of interest and the lowest intercorrelation with the other secondary variables (Table 3.). Additionally, the attribute selection was completed using a stepwise regression that results in a ranked list of variables in order of increasing evaluation power. A 10-fold cross validation was used to ensure the data was not overfit by the model.  Table 3. Ranked results of a filter-based attribute selection over the full set of data used from over Turrialba volcano.

Variable Rank
Soil Temperature (°C) 1 Air Temperature (°C) 2 Relative Humidity 3.2 Soil Moisture (VWC) 3.9 Pressure (mBar) 4.9 Elevation (m) 6 Each of 12 regression models were developed to analyze the interconnection between CO2 gas flux and external meteorological variables. Based on the results of the attribute selection that ranked the secondary variables, each model was run six times, each time adding a variable in order of decreasing rank. The best model overall was selected as that which has the highest correlation coefficient and the lowest mean absolute error. Again, a 10-fold cross validation was used in every case to avoid overfitting the model. The best model over these data was a linear regression model using all six variables. The coefficient of correlation is 0.25, and the mean absolute error is 35.67, indicating there is limited predicting power of this model. Full results are displayed in Appendix A-1.
To test the best model trained over the full dataset the 27 measurements at the control point each day of measurement are fed into the same model within the WEKA interface. Using the same model for both data sets avoids any error introduced by normalization that may be created if considering each data set independently, especially as the difference in the number of measurements is so great. The results of this model over the testing data show a coefficient of correlation equaling 0.17 and a mean absolute error of 1.04. This indicates that the model fits the data well but does not provide much significant predicting power. These results are in disagreement with previous studies showing influence on CO2 flux from barometric pressure and humidity, in particular, though these were exacerbated by changes between seasons and diurnally (Chiodini et al., 1998;Hinkle et al., 1994).
We hypothesize that the periods over which we are measuring are relatively short, and so do not adequately capture those variations in meteorological conditions that are great enough to significantly alter the amount of fCO2. Data over long time-frames, particularly over multiple seasons, would be required to accurately model the influence on fCO2 measurements. Measurement of the variation related to the diurnal effect on fCO2 would require analysis outside the scope of this study. It is possible that in low flux areas, such as the control point, small changes in fCO2 are related to natural variability in the CO2 diffusion gradient rather than from a significant influence from measured external variables. There is also potential for small variations in soil conditions, or a combination of meteorological conditions, to create variability in fCO2 between measurement times. If there is an influence that we are not aware of we stipulate the error would be no greater than 10%, and relative variability in the data would be preserved. For this study, the effect of any meteorological influence not measured are treated as negligible under the assumption that the influence would not be so great as to change the relative population a flux measurement falls within.

δ 13 C Isotope Populations
Isotopic concentrations, in parts per mille (per thousand), were collected over each area during fCO2 measurements. Results of the 53 samples are displayed in Figure 6. with a range of values calculated by Malowany et al., (2017) representing the volcanic end member based on plume samples and fumarole samples taken from Turrialba volcano in 2012. Values of δ 13 C fall between -30‰ and 0‰ and are divided by the geographical area they were taken from. The highest δ 13 C values collected were from the summit of Turrialba (maximum -1.19‰) and the lowest were collected from the Santa Cruz area (Fig. 6) (minimum -25.91‰). The average error is 1.8‰, with higher error derived from more diluted samples (full data displayed in supplementary material). Methane gas is a useful tool toward determining carbon flux origins, even as it is impossible to clearly define the origin of methane itself. There are four possible origins of methane emissions in our collected samples: 1) introduction as atmospheric contamination during sample collection; 2) air-saturated ground or meteoric water; 3) production by the activity of microbes or other organisms in the surrounding soil, or 4) introduction by hydrothermal fluid and vapor interactions (Cerling, 1984;Etiope et al., 2006;Giggenbach, 1996). The ratio between CH4/CO2 further constrains the origin of gas flux, even as the origin of methane is poorly constrained. The ratio between CH4 and CO2 from each sample were calculated from CH4 and CO2 concentrations (in ppm). Values range from 2.05E10 -5 at the summit of Turrialba to 0.0044 on the flanks of Irazú, as shown in Figure 7. The average error is 0.001 for these measurements, with higher error from more diluted samples. Isotopic values do not show any significant correlation across the dataset either with CO2 flux values or the commonly used inverse of CO2 flux values. Because of small scale spatial heterogeneity within CO2 flux it is possible that though CO2 isotopes were collected within 1 m of the accumulation chamber this may be a great enough distance to cause a significant change in diffuse degassing. It may be more likely, however, that the origin of CO2 gas is not directly dependent on the intensity of CO2 flux present. Assuming this is the case, proportions of each origin of CO2 gas (indicated isotopically) would remain relatively the same as when CO2 from the melt mixed with other end members at the source or along ascent. This would explain the lack of any correlation between the amount of degassing (CO2 flux) and isotopic composition (δ 13 C). The only exception where we do see a significant relationship is between δ 13 C and 1/fCO2 over the flanks of Irazú volcano (Fig.  2). Here, there is an R 2 value of 0.61, indicating a significant correlation. There were no values of high flux in this area, the land cover is predominantly agricultural, and isotopic values are generally depleted. Therefore, we suggest that the main driver of fCO2 is from a biogenic source that has a relatively low contribution from a magmatic source. There is no generalized correlation between isotopic values and ground temperature or elevation. The lack of correlation between δ 13 C and ground temperature is easily explained by a stronger influence from the sun's heat than any subsurface magma, especially where volcanically-derived fCO2 is not present. Ground temperature does not correlate with fCO2 until a high flux and high temperature threshold is reached. The actual value will vary between studies and volcanic systems, but the same phenomena is visible across many studies (Lin et al., 2019, Werner et al., 2014. Our primary interest is in diffuse degassing, meaning areas with temperatures and CO2 flux concentrations above this threshold would likely be outside the scope of this study. Topographic controls on degassing are expected in the study area. For example, areas of high intensity fCO2 along the Ariete fault are concentrated at the base of a fault scarp, which also acts as an intermittent stream. However, to explore any correlation numerically would necessitate use of relative elevation (or relief) rather than absolute elevation as the study area spans more than 2000 m in elevation.

Gas Distribution 2013 to Present
Previous The hydrothermal reservoir present between the shallow magma reservoir and the summit area of Turrialba is alternately boiled-off and replenished, and magmatic gas does not always travel through the reservoir but rather can bypass it. We postulate that measurements of CO2 flux at the summit likely represent contributions from both the shallow and deep magma reservoirs, following the conclusions of Camarda et al. (2012). Any changes to summit degassing between 2013 and the present may reflect changes in the geometry of the magma chamber, the permeable pathways gas may travel along toward the surface and/or changes within the magmatic or hydrothermal system. Our measurements across the broader edifice provide insight into the extent of the magmatic system, representing both deep and shallow magmatic CO2 sources (Camarda et al., 2012).
We completed two days of measurements, resulting in 82 measurement points, over representative cross sections of the summit area to compensate for limited time and access. Three areas were identified as being geologically unique based on historic activity and average intensity and, therefore, were analyzed independently. These are a fumarolic field to the southwest of the west active crater, the central crater, and the eastern crater (Fig. 8).
The highest intensity degassing is observed in the west summit crater, recognized as the source of a plume intermittently visible above the volcanic summit. Gas contributions from this area are not considered in this study as emissions are above the limit of the sensor and physically inaccessible. We use direct comparisons between 2013 and 2021 data collection points within the maximum horizontal accuracy of a standard GPS unit (~10 m) to examine potential changes in the summit gas emissions through time (Table 4.). Turrialba has shown significant eruptive activity in the eight years between the Epiard et al. (2017) and present measurement campaigns. Generally, degassing seems to be continuing along a multi-year trend toward the West-Southwest, leaving the Eastern crater devoid of observable degassing. The average flux over the Eastern crater was low to begin with (3.6 g /m 2 /day), suggesting the average decrease in flux since 2013 of 1.7 g/m 2 /day to be a significant change. Degassing anomalies at the summit are concentrated along the SW-NE trending lineaments that are part of the Turrialba summit graben (Fig. 9) (Epiard et al., 2017;Montero et al., 2011). Summit structures have been defined as a series of normal faults connecting the Elia fault, which runs along the northern edge of the collapse structure encompassing the NE edge of the Turrialba edifice to the Ariete fault along the southern flank (Fig. 2) (Calvo et al., 2019;Montero et al., 2011). Both the Elia and Ariete faults trend SW-NE, while the connecting faults trend slightly more to the north. Extension and a pull-apart structure opened a pathway connecting mantle-derived melt with the surface, allowing initial formation of Turrialba and Irazú (Montero et al., 2011). Normal motion along a structure has been linked to degassing where rifting is present (Jolie et al., 2019) and over active volcanoes (De Gregorio et al., 2002). Historical motion along the summit normal and strike-slip faults with a normal component suggest the potential for a similar process here, leading to summit degassing features that follow these structures. Since the most recent measurement campaign in 2013, degassing continues to be concentrated along one lineament with fumarolic fields aligned along the northeast and southwest sides of the active west crater. Each unique area of analysis shows a decrease in diffuse degassing, with the greatest change being within the central crater (an average of -1984.5 g m -2 day -1 ). Decreasing diffusive gas flux at the summit overall is likely due to: 1) degassing becoming more channelized, concentrating at the active crater and along this lineament (and in associated fumarolic fields) or, 2) a thick blanket of ash produced by recent eruptions is obstructing gas flow. It is possible that a combination of both is influencing the distribution of gas flow; though more analysis would be required to determine which may be the dominant cause with greater confidence. Additionally, areas along this lineament show evidence of long-term alteration that may lead to greater instability and greater likelihood of collapse in the future. For example, in October 2014 a highly altered area on the western edge of the western crater collapsed to form the larger crater still present today (Epiard et al., 2017).

Carbon Isotopic Signatures
The signature of carbon isotopes collected over the study area shows relative contributions from each of three carbon source end members. Mixing between volcanic, biogenic, and atmospheric sources define the range of δ 13 C and CH4/CO2 values, all of which are utilized in this study (Fig. 10). Carbon 13 isotopes are often used as an indicator of volcanic activity as they are easily distinguished from carbon 12, which is the primary component of atmospheric carbon species (Chiodini et al., 2008;Malowany et al., 2017;Vaselli et al., 2009;and others). One primary source of 13 C globally is from volcanic activity, where atmospheric carbon is considerably depleted in 13 C by comparison (NOAA). Because there is almost no δ 13 C respired biotically, any 13 C present is likely due to volcanic activity or is residual in the atmosphere where there is no significant industrial contribution (Bogue et al., 2019). Atmospheric δ 13 C values fall between biogenic and volcanic δ 13 C, meaning it can be difficult to identify the atmospheric component alone based only on δ 13 C values. Methane, and further, the ratio of methane to carbon dioxide in a sample allow further differentiation between atmospheric δ 13 C, biogenic, and volcanic δ 13 C endmembers. Methane concentrations in the atmosphere are higher on average than those present where volcanic or biogenic emissions are dominant. Therefore, samples with high CH4/CO2 ratios are expected to have a higher atmospheric component and be "pulled" to the right when plotting δ 13 C and CH4/CO2 (Fig. 10). Atmospheric δ 13 C values are well studied and frequently measured due to their importance for carbon cycle monitoring and tracking of industrial CO2 output. The isotopic ratio of the atmosphere varies between biomes, seasonally, and diurnally on small spatial scales (NOAA). However, this study uses a generalized atmospheric δ 13 C value of -8‰ based on two assumptions. First, because all samples were taken within the same timeframe each day and within the same season any diurnal or seasonal effects would be avoided. Second, Malowany et al., (2017) recorded ambient δ 13 C values of -8.3‰ to -10.2‰ in 2014 at Turrialba volcano, which are reasonably close to the average value used in this study. As the isotope samples collected and analyzed in this study were collected from the soil, it is not possible to measure an uncontaminated "true" atmospheric end member. Instead, an average value of atmospheric δ 13 C collected by the National Oceanic and Atmospheric Administration (NOAA) is used even as the true atmospheric end member is a range. The same approach is used when finding a representative value of CH4/CO2. This was calculated on average for the atmosphere from the average concentration of CO2 in the atmosphere in 2019 (409.8 ± 0.1 ppm) and the same for methane in 2018 (1.85 ppm). Any confounding influence from industrial CO2 is considered negligible as the trade winds through this part of Costa Rica are consistently directed west southwest, from the volcano towards the Central Valley, which contains the only major urban center near Turrialba (Bogue et al., 2019).
Similarly, the biogenic isotopic end member is highly variable depending on ecological conditions, climate, and time of day. Soils have a higher partial pressure than the atmosphere due to the respiration (or release of CO2) of roots and other microbes in the soil (Cerling, 1984). Much like volcanic CO2, biotic CO2 is transported to the surface; though usually at a lower rate. The isotopic composition of organic soil in equatorial climates typically consists of respired CO2, where δ 13 C concentrations usually fall below ~-20‰ (Cerling, 1984;Malowany et al., 2017). Our data falls within the realm of known values, with a biological δ 13 C end member of 26‰ to -24‰. While plants are known to produce methane gas, estimates of concentration remain relatively low on a global scale (Houweling et al., 2006). Additionally, where methane gas remains in the soil for long periods and is not facilitated to the surface by high levels of gas flux (such as those form hydrothermal or magmatic systems) it is oxidized by the local ecosystem and not evident at the surface (Tassi et al., 2013). Finally, we expect the concentrations of CO2 respired by plant photosynthesis to overwhelm concentrations of methane, resulting in a low CH4/CO2 ratio where a biogenic source is present.
The magmatic end member composition of δ 13 C is variable among individual volcanic systems and through eruptive cycles even within the same system over time Malowany et al., 2017;Vivieros et al., 2020). The composition of surface emanations is influenced by both the source magma and overlying hydrothermal system as carbon species rise to the level of volcanic activity. Higher, more enriched values indicate a greater magmatic contribution (Chiodini et al., 2008). Typically, enriched values are representative of melt-derived CO2 and have δ 13 C values above -5‰. The cyclic behavior between hydrothermal/phreatic and magmatic/phreatomagmatic eruptive activity at Turrialba allows us to utilize past studies over periods of well constrained activity to identify the volcanic source end member composition from a time when activity was similar to now (DeMoor et al., 2016;Vaselli et al., 2009). Vaselli et al. (2009) reports δ 13 C, CO2 concentrations, and methane concentrations from two summit fumaroles, one high temperature vent, and one fumarole within the Ariete fault, respectively, from 1999 to 2008 (Fig. 1). Their results provide the basis for determining end members within the volcanic system of Turrialba as they collected high-temperature samples where we would expect the least contamination from biological or atmospheric carbon .
The primary driver of methane production in areas of intense volcanic activity is through generation in the hydrothermal reservoir (Giggenbach et al., 1996;Vaselli et al., 2009). Therefore, we use data from Vaselli et al (2009) over well-constrained active periods to identify a reasonable end member for CH4/CO2 in magmatic gas at Turrialba. There is an observable downward trend in CH4/CO2 at the west and central craters from 1999 until around 2002, when they determined the system transitioned from a dominantly hydrothermal to hydrothermal magmatic type of volcanic activity. This trend continued between 2002 and 2005, when the hydrothermal system appears to have largely dried out due to the presence of high temperature magma at shallow levels, therefore, shifting the gas composition towards a dominantly magmatic signature (Malowany et al., 2017;Vaselli et al., 2009). Periods of hydrothermal/magmatic or magmatic-dominated activity have characteristically low CH4/CO2 values as the hydrothermal signal is increasingly diluted by magmatic CO2 input (Giggenbach, 1996). The west high temperature (HT) vent (shown in Figure 1.) shows the same process on an accelerated time scale (Fig. 11). There is a decrease in CH4/CO2 as the vent, or high-temperature fumarole, increased rapidly in temperature from 2007 to 2008, which was interpreted as a shallow hydrothermal aquifer boiling off (Vaselli et al., 2009). March 2008 represents the clearest evidence of a "pure" magmatic end member within a period of measurement, and assuming the gas source of the volcanic system has remained stable over time, it is likely a good representation of a "pure" magmatic end member today. A comparison between 2008 and 2021 shows a lower CH4/CO2 ratio in 2008 than any sample at present (Fig. 12). Additionally, δ 13 C is higher at the west HT vent than at any location measured in 2021. Contrary to Vaselli et al. (2009) and other past studies, not all our samples were collected from intensely active areas (i.e., fumaroles), but instead focused on the soil slightly outside of those features. This may introduce some atmospheric contamination and cause δ 13 C values to appear more depleted than they may be within the  . Nevertheless, isotopic measurements at the summit still suggest a connection with a deep magmatic reservoir. In particular, degassing is concentrated along the SW-NE trending lineament running through the active crater and two fumarolic fields, where the most enriched (volcanic) δ 13 C values and highest diffuse CO2 fluxes are found (Fig. 2). Diffuse soil fCO2 measurements and δ 13 C values are consistent with the conclusion of Epiard et al., (2017) suggesting that there is likely gas flow from depth concentrated along this structure, possibly facilitating further concentration of degassing toward a centralized location.

Hydrothermal Component: The Ariete Fault
When hydrothermal activity is present at the summit of Turrialba, which it is periodically, it introduces slightly depleted δ 13 C values and slightly more elevated CH4/CO2 values relative to those that are related to a more magma dominant system. Previous studies have used the Ariete fault fumaroles as an indicator of increased hydrothermal activity relative to the summit (Malowany et al., 2017;Vaselli et al., 2009). Fumarole δ 13 C samples, taken from an area known for consistent fumarolic activity in the Ariete fault, retain a consistent d13C signature over many years even as summit δ 13 C values are more variable (Table 5.). Table 5. δ 13 C values collected from the Ariete fault over multiple independent measurement campaigns.

Date δ 13 C Reference
Mar-08 -2.8 Vaselli et al., 2009Sep-13 -2.96 Epiard et al., 2017Jan-14 -3.33 Malowany et al., 2017 This study The Ariete fault forms a scarp trending SW-NE, significant enough to be visible by satellite, with an intermittent stream at the bottom of a steep valley (40-50 m wide) at the base of the scarp (Fig. 2). Measurements within this valley show high fCO2 values leading toward advective degassing within the previously mentioned fumaroles. Fumarolic activity culminates at a dry waterfall when moving up the drainage, toward the summit. The zone of elevated gas flux within the drainage, including high soil flux and fumarolic activity, is bounded on either end by steep dry waterfalls. This dry waterfall possibly indicates a tectonic nickpoint, and if so, may also indicate a cross-cutting geologic structure impacting degassing pathways. A similar structure is visible near the southern edge of significant degassing along the Ariete fault. Soil gas isotopic samples taken perpendicular to the Ariete fault show a gradient from enriched δ 13 C values toward a more negative, depleted isotopic signature with increased distance from the fault (R 2 = 0.84) (Fig. 13).

More volcanic
More atmospheric would expect from magmatically-dominated gas (Fig. 14). Giggenbach (1996) determined that magmatic gas would have a very small methane component, if any at all. The presence of higher methane concentrations in the Ariete fault, therefore, are interpreted to indicate the presence of a shallow hydrothermal aquifer contributing to degassing here that is likely independent from activity at the summit.

Hazard Implications: Turrialba Western Slope
On the western slope of Turrialba lies La Silvia Quemada, the burned forest, named as such because it is the site of mass vegetation death, caused by strongly elevated soil and plume volcanic gas levels combined with ash fallout since 2004 (Fig. 2) (Tortini et al., 2017). Additionally, there are areas of diffuse degassing so intense that burned vegetation appeared toward the beginning of activity high on the NW flank, but to our knowledge there have been no diffuse soil degassing studies covering this area (Global Volcanism Program, 2013). We record a lack of significant fCO2 throughout this area, with the lowest values of any area studied (Table 1.). Isotopic signatures and CH4/CO2 values fall between all three end members (Fig. 10), making any interpretation of activity difficult. A thick layer of ash covers the entire area, with the potential to reduce gas permeability and subsequently measurable flux. Fumarolic activity was observed, however, just below the outer cone of the active crater near the origin of a 2016 lahar. Field conditions made measurements in this area impossible, limiting our analysis to a visual recognition of advective degassing. Epiard et al., (2017) have hypothesized a weakening of the upper outer cone, consistent with this observation. Future studies should consider weakening from chemical alteration as a component of hazard monitoring at Turrialba, especially in the scenario where gas is unable to pass through a layer of semi-impermeable ash and instead pools beneath it or is channelized through another structure. If there is a semiimpermeable layer between the subsurface and diffuse degassing measurement equipment it is important for hazard prediction and monitoring to find an alternative method when determining the extent of Turrialba's magmatic activity to the west.

Activity between Turrialba and Irazú
In exploring the potential connection between Turrialba and Irazú, or if there is an independent presence of magma at depth, we are interested in any diffuse degassing anomalies in the depression between the two volcanoes, the San Gerardo area, or across the flanks of Irazú and Turrialba (the Santa Cruz area) (Fig. 2). The Irazú and Santa Cruz areas have characteristically low flux values, consistent with a lack of known faults or features in either area. There are no documented faults or fractures within the San Gerardo either, and there have not been any diffuse degassing studies in this depression to date. However, there is an actively degassing hot spring (the Santa Teresa spring) and additional warm-and cold-CO2 springs within the banks of a river that were included in this study. There were anomalously high fCO2 values associated with all three springs investigated (a maximum flux at Santa Teresa of 540.87g m -2 day -1 , and maximum fluxes of 232.7 g m -2 day -1 , 29.5 g m -2 day -1 , and 34.5 g m -2 day -1 of three additional springs, respectively).
Isotopic 13 C values from San Gerardo fall along a well-defined linear mixing line (R 2 =0.8) between biogenic and atmospheric end members, with samples within this trend predominantly being collected from the Irazú, Santa Cruz, and San Gerardo areas (Fig. 15). Initially, this would suggest that these three areas are unlikely to have been influenced by volcanic gas input and instead reflect atmospheric contamination during sampling, which increases with higher dilution of the sample. However, the intensity of CO2 degassing from the San Gerardo area is inconsistent with a biogenic or atmospheric source. The difference between the isotopic signature and intensity fCO2 from what we would expect may be related to a greater than normal distance between fCO2 measurements and isotope sample collections, though unlikely. Samples collected in-situ from the Santa Teresa spring by OVSICORI in 2017 reflect conditions impossible to measure during the 2021 sampling campaign. This spring has moderate temperatures (~40-50°C) and is characterized by mineral deposits typical of hydrothermal activity, such as hematite and goethite, and red precipitates in the water (location displayed in Fig. 1). Rouwet et al. (2021) obtained chemistry of the water and precipitates from the Santa Teresa spring in 2007, 2009, and 2012. The δ 13 C and CH4/CO2 values from 2017 fall within a reasonable range of other known hydrothermal activity at Turrialba, such as the Ariete fault (Fig. 14). This is consistent with spring chemistry, indicating steam-heating and near-surface vapor-and-liquid interactions (Rouwet et al., 2021).
There are multiple possible origins for the Santa Teresa hot spring and other CO2 springs nearby. Their low elevation relative to the volcanic centers of both Irazú and Turrialba and their proximity to an active hydrologic system suggest the possibility that volcanic CO2 is being captured farther up in the volcanic system and washed down before resurfacing at the base of this depression. Alternately, such high levels of fCO2 suggest the presence of previously unmapped melt at depth degassing to the surface through a similarly unmapped DDS. This could indicate from the presence of an independent body of magma, or a magma plumbing system that connects those underlying and directly feeding the Irazú and Turrialba volcanoes. In order to fully evaluate the area, or to determine if there is in fact magma present at depth below the depression, we would need to expand the area surveyed by diffuse degassing measurements and isotopic gas sampling. Collections of water chemistry data, geophysical imaging, and measurements of deformation would also greatly increase certainty when modeling this area and the volcanic systems of both volcanoes.

Conclusion
Our investigation of soil diffuse CO2 degassing has shown the following:

I.
Isotopic CO2 values at the summit of Turrialba reflect those made in past studies (Malowany et al., 2017;Vaselli et al., 2009), and support a deep degassing connection between the lineament of a summit graben and a deep magmatic reservoir (Epiard et al., 2017). Concentrations of δ 13 C remain relatively consistent over the years where data are present. II.
The 2014 collapse of the eastern wall into the western crater combined with visual observations of degassing under the edge of the upper outer cone and continued high-intensity degassing in the fumarolic field south of the active crater support inferences that chemical alteration is continuing to destabilize the summit area of Turrialba. Though not recorded in this study, other tectonic structures interpreted as part of the summit graben have potential to begin following the same behavior, with serious hazard implications. III.
Focusing on large-scale diffuse degassing patterns, rather than tightly spaced measurements over a small feature, require a unique set of considerations. The commonly used method of measuring in a grid pattern for applications of a statistical interpolation map is not possible on this scale. Instead, it is more time and cost efficient to concentrate measurements around areas of interest (i.e. mapped faults and features) throughout the study area, with "background" measurements in-between. Maps and past studies are essential to this process, along with the knowledge of those who live and work in the area. Their contributions are crucial when discussing locations and timelines of volcanic activity. IV.
The ratio of methane to carbon dioxide has proven a useful tool toward determining the source mechanism and contribution of diffuse soil CO2 flux between biogenic, atmospheric, and magmatic end members, with potential to be applied regarding differentiation between magmatic and hydrothermal-magmatic degassing. Because the same analyzer can be used both on CO2 and CH4, the ratio between CH4 and CO2 has potential act as a cost-effective and convenient alternative to more difficult to sample volcanic gases often used to track gas origin, such as helium or radon. We recommend further analysis include such considerations and test this conclusion across other volcanic systems.
This is the first study at Turrialba volcano where: a) soil gas emissions are considered significantly beyond the crater vent area, helping to distinguish between vent fluxes, fault degassing, and background degassing, and b) measurements were made to investigate a connection between the Turrialba and Irazú volcanic systems across the valley that separates them.

Future Directions
Implications from diffuse, elevated CO2 concentrations emitted by volcanic systems in tropical climates extend from changing physiological processes of tropical trees to allowing the monitoring of subsurface melt as it pertains to eruptive hazards and connectivity between volcanic systems. Future directions would include both of these avenues. Wood carbon isotopes have proven to reflect uptake of volcanic carbon as the tree grows and matures (Bogue et al., 2019). Further constraints on the path of fCO2 after it is released from the ground and on how much volcanic carbon is absorbed by a tree in a given area allow delineation of CO2 diffusion and provide insight about terrestrial carbon cycling. The area between Turrialba and Irazú is of particular interest. The presence of elevated fCO2 from springs throughout the depression imply there is a volcanic source feeding them.
Further measurements in this area are essential for understanding if there is melt at depth which could cause future volcanic hazards, and toward determining if Turrialba and Irazú have any connection between their magma reservoirs. Diffuse degassing campaigns are time consuming, and the resolution of the resultant data is often too low to accurately interpolate carbon flux over large areas. Cawse-Nicholson et al., (2019) and others have proven remote sensing indices to show changes in CO2 flux concentrations over Mammoth Mountain, USA by using the boreal forest as a predictor. Further analysis is required to apply these same techniques in a tropical forest, but the result would be an efficient tool for measuring volcanic CO2 over large temporal and spatial scales.

A Error Regression Analysis
Only 395 measurements points were possible to use in the multivariate analysis as they all contained measurements of all variables of interest (results in Table A-1). Points with missing meteorological measurements were not considered in the error analysis. The most common cause of a missing variable is hard ground which is impenetrable to soil temperature and moisture probes. 27 values were recorded at a control point on the southern flanks.  1955 31.1947 0.1955 31.1947 0.1955 31.1947 0.1955 31.1947 0.1955 31.1947 0.1955 31.1947