Plant responses to volcanically elevated CO₂ in two Costa Rican forests

We explore the use of active volcanoes to determine the short- and long-term effects of elevated CO₂ on tropical trees. Active volcanoes continuously but variably emit CO₂ through diffuse emissions on their flanks, exposing the overlying ecosystems to elevated levels of atmospheric CO₂. We found tight correlations (r²=0.86 and r²=0.74) between wood stable carbon isotopic composition and co-located volcanogenic CO₂ emissions for two of three investigated species (Oreopanax xalapensis and Buddleja nitida), which documents the long-term photosynthetic incorporation of isotopically heavy volcanogenic carbon into wood biomass. Measurements of leaf fluorescence and chlorophyll concentration suggest that volcanic CO₂ also has measurable short-term functional impacts on select species of tropical trees. Our findings indicate significant potential for future studies to utilize ecosystems located on active volcanoes as natural experiments to examine the ecological impacts of elevated atmospheric CO₂ in the tropics and elsewhere. Results also point the way toward a possible future utilization of ecosystems exposed to volcanically elevated CO₂ to detect changes in deep volcanic degassing by using selected species of trees as sensors.

1 Introduction

Tropical forests represent about 40 % of terrestrial net primary productivity (NPP) worldwide, store 25 % of biomass carbon, and may contain 50 % of all species on Earth, but the projected future responses of tropical plants to globally rising levels of CO₂ are poorly understood (Leigh et al., 2004; Townsend et al., 2011). The largest source of uncertainty comes from a lack of understanding of long-term CO₂ fertilization effects in the tropics (Cox et al., 2013). Reducing this uncertainty would significantly improve Earth system models, advances which would help better constrain projections in future climate models (Cox et al., 2013; Friedlingstein et al., 2013). Ongoing debate surrounds the question of how much more atmospheric CO₂ tropical ecosystems can absorb; this is known as the “CO₂ fertilization effect” (Gregory et al., 2009; Kauwe et al., 2016; Keeling, 1973; Schimel et al., 2015).

Free Air CO₂ Enrichment (FACE) experiments have been conducted to probe this question, but none have been conducted in tropical ecosystems (e.g., Ainsworth and Long, 2005; Norby et al., 2016). Some studies have used CO₂-emitting natural springs to study plant responses to elevated CO₂, but these have been limited in scope due to the small spatial areas around springs that experience elevated CO₂ (Hattenschwiler et al., 1997; Körner and Miglietta, 1994; Paoletti et al., 2007; Saurer et al., 2003). These studies have suffered from several confounding influences, including other gas species that accompany CO₂ emissions at these springs, human disturbances, and difficulty with finding appropriate control locations. Additionally, none have been conducted in the tropics (Pinkard et al., 2010). A series of studies in Yellowstone National Park (USA) used its widespread volcanic hydrothermal CO₂ emissions for the same purpose, though it is not in the tropics (Sharma and Williams, 2009; Tercek et al., 2008). Yellowstone was particularly suitable for this type of study due to its protected status as a national park and because the large areas of CO₂ emissions made control points more available (Sharma and Williams, 2009; Tercek et al., 2008). These studies reported changes in rubisco, an enzyme central to CO₂ fixation, and sugar production in leaves, similar to results from FACE experiments, suggesting that volcanically influenced areas like Yellowstone have untapped potential for studying the long-term effects of elevated CO₂ on plants.

Tropical ecosystems on the vegetated flanks of active volcanoes offer large and diverse ecosystems that could make this type of study viable. Well over 200 active volcanoes are in the tropics (Global Volcanism Program, 2013), and many of these volcanoes are heavily forested. However, fewer of these tropical volcanic forests have sufficient legal protection for being a source of long-term information, and the effects of diffuse volcanic flank gas emissions on the overlying ecosystems remain largely unknown. Most previous studies focused on extreme conditions, such as tree kill areas associated with extraordinarily high CO₂ emissions at Mammoth Mountain, CA (USA; Biondi and Fessenden, 1999; Farrar et al., 1995; Sorey et al., 1998). However, the non-lethal effects of cold volcanic CO₂ emissions – away from the peak emission zones, but still in the theorized fertilization window – have received little attention and could offer a new approach for studying the effects of elevated CO₂ on ecosystems (Cawse-Nicholson et al., 2018; Vodnik et al., 2018). The broad flanks of active volcanoes experience diffuse emissions of excess CO₂ because the underlying active magma bodies continuously release gas, dominated by CO₂ transported to the surface along fault lines (Chiodini et al., 1998; Dietrich et al., 2016; Farrar et al., 1995). This process has frequently been studied to understand the dynamics of active magma chambers and to assess potential volcanic hazards (Chiodini et al., 1998; Sorey et al., 1998). These emissions are released through faults and fractures on the flanks of the volcano (Burton et al., 2013; Pérez et al., 2011; Williams-Jones et al., 2000; see Supplement Fig. S1). Volcanic flanks through which these gases emanate are broad, typically covering 50–200 km², and often support well-developed, healthy ecosystems. Some of these faults tap into shallow acid hydrothermal aquifers, but by the time these gases reach the surface of most forested volcanoes, soluble and reactive volcanic gas species (e.g., SO₂, HF, HCl, and H₂S) have been scrubbed out in the deep subsurface, leading to a diffusely emanated gas mix of predominantly CO₂ with minor amounts of hydrogen, helium, and water vapor reaching the surface (Symonds et al., 2001).

Trees in these locations are continuously exposed to somewhat variably elevated concentrations of CO₂ (eCO₂), although the specific effects of this eCO₂ on the trees are not well understood. Volcanic CO₂ has no ¹⁴C and a δ¹³C signature typically ranging from around −7 ‰ to −1 ‰, which is distinct from typical vegetation and noticeably enriched in ¹³C compared to typical atmospheric values (Mason et al., 2017). If trees incorporate volcanic CO₂, then the stable carbon isotopic composition of wood may document the long-term, possibly variable influence of volcanic CO₂ during the tree's growth. With this tracer available, volcanic ecosystems could become a valuable natural laboratory for studying the long-term effects of elevated CO₂ on ecosystems, especially in understudied regions like the tropics. Several studies have found correlations between variations in volcanic CO₂ flux and plant ¹⁴C records at Mammoth Mountain, Yellowstone, and Naples, which agreed well with previous observations at these well-studied sites (Cook et al., 2001; Evans et al., 2010; Lefevre et al., 2017; Lewicki et al., 2014). The Mammoth Mountain and Yellowstone studies linked seismic swarms and accompanying increases in CO₂ flux to decreases in ¹⁴C content in tree rings in one or two trees, demonstrating the method's utility for uncovering yearly scale variations in volcanic CO₂ fluxes (Cook et al., 2001; Evans et al., 2010; Lewicki et al., 2014). The Naples study instead focused on using ¹⁴C in grasses as short-term (2 to 6 months) monitors of volcanic CO₂ flux, which is useful for volcanic monitoring due to the time-integrated signal they provide (Lefevre et al., 2017). A study of plants growing at Furnas volcano found very strong (r²>0.85) correlations between depletions in ¹⁴C and enrichments in ¹³C from volcanic CO₂ in three species of plants, although this study also had a relatively limited (five samples per species) dataset (Pasquier-Cardin et al., 1999). The previously mentioned Naples study also found some correlation between ¹³C and ¹⁴C, although it was not as strong as the study in Furnas (Lefevre et al., 2017; Pasquier-Cardin et al., 1999). Additionally, short-term effects of eCO₂ might be revealed by plant functional measurements at the leaf scale, where the additional CO₂ could increase carbon uptake in photosynthesis. A series of studies at Mt. Etna in Italy and Mt. Nyiragongo in the Democratic Republic of the Congo found linear anomalies in the NDVI (normalized difference vegetation index), a measure of vegetation greenness (Houlié et al., 2006). One to two years after the appearance of the NDVI anomalies, flank eruptions occurred directly along the line of the anomaly, indicating a plant response to the volcano's pre-eruptive state which may be due to increased CO₂ emissions in the buildup to the eruption (Houlié et al., 2006). A follow-up study found that the trees on Mt. Etna were relatively insensitive to changes in temperature and water availability, strengthening the case that volcanic influence was indeed responsible for the NDVI anomaly (Seiler et al., 2017).

Here we provide preliminary results on the short- and long-term non-lethal impacts of diffuse volcanic CO₂ emissions on three species of tropical trees on the flanks of two active volcanoes in Costa Rica. We also explore the viability in studying volcanically influenced ecosystems to better understand potential future responses to elevated CO₂ and suggest adjustments to our approach that will benefit future, similarly motivated studies.

Figure 1 Overview of measurement locations in two old-growth forests on the upper flanks of two active volcanoes in Costa Rica, Turrialba and Irazú. Distribution of mean soil CO₂ flux across north flank of Irazú (a) and south flank of Turrialba (b). Colors of dots correspond to flux populations (see Fig. 3).

2 Methods

2.1 Investigated locations and sampling strategy

Irazú and Turrialba are two active volcanoes located ∼25 and 35 km east of San José, Costa Rica (Fig. 1). These two volcanoes are divided by a large erosional basin. The forested portions of the two volcanoes cover approximately 315 km². The vast majority of the northern flanks of Irazú and Turrialba are covered in legally protected dense old-growth forest, while the southern flanks are dominated by pastureland and agriculture. Turrialba rises 3300 m above its base and has been active for at least 75 000 years, with mostly fumarolic activity since its last major eruption in 1866 (Alvarado et al., 2006). It has experienced renewed activity beginning in 2010, and its current activity is primarily characterized by a near-constant volcanic degassing plume, episodic minor ash emissions, and fumarolic discharges at two of the summit craters, as well as significant diffuse and fumarolic gas emissions across its flanks, focused along fault systems (Martini et al., 2010). Turrialba's CO₂ emissions in areas proximal to the crater were calculated at 113±46 tons d⁻¹ (Epiard et al., 2017). The Falla Ariete (Ariete fault), a major regional fault, runs northeast–southwest through the southern part of Turrialba's central edifice and is one of the largest areas of diffuse CO₂ emissions on Turrialba (Epiard et al., 2017; Rizzo et al., 2016). Atmospheric CO₂ has an average δ¹³C value of −9.2 ‰ at Turrialba, and the volcanic CO₂ released at the Ariete fault has significantly heavier δ¹³C values clustered around −3.4 ‰ (Malowany et al., 2017).

Irazú has been active for at least 3000 years and had minor phreatomagmatic eruptions in 1963 and a single hydrothermal eruption in 1994. Currently, Irazú's activity primarily consists of shallow seismic swarms, fumarolic crater gas emissions, small volcanic landslides, and minor gas emissions on its northern forested flank (Alvarado et al., 2006; Barquero et al., 1995). Diffuse cold flank emissions of volcanic CO₂ represent the vast majority of gas discharge from Irazú, as the main crater releases 3.8 t d⁻¹ of CO₂ and a small area on the north flank alone releases 15 t d⁻¹ (Epiard et al., 2017). Between the two volcanoes, a major erosional depression is partially occupied by extensive dairy farms and is somewhat less forested than their flanks.

In this study, we focused on accessible areas between 2000 and 3300 m on both volcanoes (Fig. 1). On Irazú, we sampled trees and CO₂ fluxes from the summit area to the north, near the approximately north–south striking Rio Sucio fault, crossing into the area dominated by dairy farms on Irazú's lower northeastern slope. Of significant importance for this type of study is that all active volcanoes on Earth continuously emit CO₂ diffusely through fractures and diffuse degassing structures on their flanks, at distances hundreds to thousands of meters away from the crater (Dietrich et al., 2016; Epiard et al., 2017), and this elevated CO₂ degassing persists continuously and consistently over decades to centuries (Burton et al., 2013; Delmelle and Stix, 2000; Nicholson, 2017). There is no inherent seasonal or meteorological variability in the source gas pressure, and there is no dependence on shallow soil or vegetation chemistry or biology (though increased soil moisture in the rainy season, wind, and atmospheric pressure can modulate gas permeability of the shallow soil; Camarda et al., 2006). The soil overlying deep-reaching fracture systems acts as a diffuser through which the volcanic gas percolates and enters the sub-canopy air. For our study sites, portions of the volcanoes with active “cold” CO₂ degassing have already been assessed and mapped previously (Epiard et al., 2017; Malowany et al., 2017).

Our sampling locations on Irazú were located along a road from the summit northward down into the low-lying area. On Turrialba, we focused on an area of known strong emissions but intact forests on the SW slope, uphill of the same erosional depression but cross-cut by the major NE–SW trending active fracture system of the Falla Ariete. We sampled three main areas of the Falla Ariete, each approximately perpendicularly transecting the degassing fault along equal altitude, at the upper Ariete fault, the lower Ariete fault, and a small basin directly east of the old Cerro Armado cinder cone on Turrialba's southwestern flank. We took a total of 51 tree samples (17 were excluded after stress screening) at irregular intervals depending on the continued availability and specimen maturity of three species present throughout the transect.

All transects are in areas experiencing measurable CO₂ enhancements from the Falla Ariete but that are not high enough in altitude to be in areas generally downwind of the prevailing crater emissions plume (Epiard et al., 2017). We avoided areas that experience ash fall, high volcanic SO₂ concentrations, and local anthropogenic CO₂ from farms or that were likely to have heavily acidified soil. Excessively high soil CO₂ concentrations can acidify soil, leading to negative impacts on ecosystems growing there (McGee and Gerlach, 1998). Because such effects reflect by-products of extreme soil CO₂ concentrations rather than direct consequences of elevated CO₂ on plants, we avoided areas with CO₂ fluxes high enough to possibly cause noticeable CO₂-induced soil acidification. Light ash fall on some days likely derived from atmospheric drift, as we were not sampling in areas downwind of the crater. The ash fall did not, in any noticeably way, affect our samples, as trees showing ash accumulation on their leaves or previous damage were the exception and avoided. Altitude, the amount of sunlight during measurements, and aspect had no consistent correlations with any of the parameters we measured.

2.2 Studied tree species

Our study focused on three tree species found commonly on Turrialba and Irazú: Buddleja nitida, Alnus acuminata, and Oreopanax xalapensis. B. nitida is a small tree with a typical stem diameter (DBH) ranging from 5 to 40 cm that grows at elevations of 2000–4000 m throughout most of Central America (Kappelle et al., 1996; Norman, 2000). The DBH of the individuals we measured ranged from 11.5 to 51.3 cm, with an average of 29.85 cm. It averages 4–15 m in height and grows primarily in early and late secondary forests (Kappelle et al., 1996; Norman, 2000). A. acuminata is a nitrogen-fixing pioneer species exotic to the tropics that can survive at elevations from 1500 to 3400 m, although it is most commonly found between 2000 and 2800 m (Weng et al., 2004). The trees we measured had DBH ranging from 14.3 to 112 cm, with an average of 57.14 cm. O. xalapensis thrives in early and late successional forests, although it can survive in primary forests as well (Kappelle et al., 1996; Quintana-Ascencio et al., 2004). It had the smallest average DBH of the three species, ranging from 6.6 to 40.9 cm, with an average of 22.71 cm.

2.3 CO₂ concentrations and soil diffuse flux measurements

Soil CO₂ flux was measured with an accumulation chamber near the base of the tree (generally within 5 m, terrain permitting) at three different points and then averaged to provide a single CO₂ flux value to compare to the ¹³C measurement of the corresponding tree sample. This technique is intended to provide a simple relative way to compare the CO₂ exposure of different trees, as a tree with high CO₂ flux near its base should experience consistently higher CO₂ concentrations than a tree with lower CO₂ flux. We also measured concentrations at ground level and 1.5–2.0 m above ground level, though these were expectedly highly variable in time and location. We analyzed CO₂ fluxes, not concentrations, because the diffuse emissions of excess volcanic CO₂ through the soil, fed from a deep magma source and dependent on the location of constant deep geological permeability, are highly invariant in time compared to under-canopy air concentrations. In contrast, instantaneous concentration measurements in the sub-canopy air are modulated by many factors, including meteorology, respiration of vegetation and animals, uptake by plants for photosynthesis, and diurnal dynamic and slope effects. An approach of instantaneous highly variable concentration measurements is thus not representative of long-term exposure. The approach of measuring the largely invariant soil-to-atmosphere volcanic CO₂ fluxes is much more representative of long-term exposure, varying mostly spatially, and the site-to-site differences are therefore more representative of the lifetime of exposure of the trees.

We used a custom-built soil flux chamber system which contained a LI-COR 840A non-dispersive infrared CO₂ sensor (LI-COR Inc., Lincoln, NE, USA) to measure soil CO₂ flux. A custom-built cylindrical accumulation chamber of defined volume was sealed to the ground and remained connected to the LI-COR sensor. The air within the accumulation chamber was continuously recirculated through the sensor, passing through a particle filter. The sensor was calibrated before deployment and performed within specifications. We recorded cell pressure and temperature; ambient pressure; air temperature; GPS location; timestamps; location description; wind speed and direction; relative humidity; and slope, aspect, and altitude as ancillary data. In typical operation, each measurement site for flux measurements was validated for leaks (visible in the live data stream display as spikes and breaks in the CO₂ concentration slope), and potential external disturbances were avoided (such as vehicle traffic, generators, or breathing animals and humans). Measurements were recorded in triplicate for at least 2 min per site. Data reduction was performed using recorded timestamps in the dataset, with conservative time margins to account for sensor response dead time, validated against consistent slope sections of increasing chamber CO₂. Fluxes were computed using ancillary pressure and temperature measurements and the geometric chamber constant (chamber volume at inserted depth, tubing volume, and sensor volume). Care was taken to not disturb the soil and overlying litter inside and adjacent to the chamber.

2.4 Leaf function measurements

Chlorophyll fluorescence measurements were conducted on leaves of all three species during the field campaign to obtain information on instantaneous plant stress using an OS30p+ fluorometer (Opti-Sciences Inc., Hudson, NH, USA). Five mature leaves from each individual tree were adapted to the dark for at least 20 min to ensure complete relaxation of the photosystems. After dark adaptation, initial minimal fluorescence was recorded (F_o) under conditions where we assume that the photosystem II (PSII) was fully reduced. Immediately following the F_o measurement, a 6000 µmol m⁻² s⁻¹ saturation pulse was delivered from an array of red LEDs at 660 nm to record maximal fluorescence emission (F_m) when the reaction centers are assumed to be fully closed. From this, the variable fluorescence was determined to be $F_{\mathrm{v}}/F_m=$ ($F_{\mathrm{m}}-F_{\mathrm{o}})/F_{\mathrm{m}}$. F_v∕F_m is a widely used chlorophyll fluorescence variable used to assess the efficiency of the PSII and, indirectly, plant stress (Baker and Oxborough, 2004). The five F_v∕F_m measurements were averaged to provide a representative value for each individual tree. Some trees had less than five measurements due to the dark adaptation clips slipping off the leaf before measurements could be taken. Ten trees had four measurements, and another six had three measurements.

The chlorophyll concentration index (CCI) was measured with a MC-100 Apogee Instruments chlorophyll concentration meter (Apogee Instruments, Inc., Logan, UT, USA). The CCI was converted to chlorophyll concentration (µmol m⁻²) with the generic formula derived by Parry et al. (2014). Depending on availability, between three and six leaves were measured for the CCI for each tree and then averaged to provide a single value for each tree. If leaves were not within reach, a branch was pulled down or individual leaves were shot down with a slingshot and collected. Photosynthetically active radiation was measured at each tree with a handheld quantum meter (Apogee Instruments, Inc., Logan, UT, USA; Supplement Table S2). Stomatal conductance to water vapor, g_s (mmol m⁻² s⁻¹), was measured between 10:00–14:00 h using a steady-state porometer (SC-1, Decagon Devices, Inc., Pullman, WA, USA), calibrated before use and read in manual mode. This leaf porometer was rated for humidity < 90 %, and humidity was sometimes above this limit during our field work. Consequently, we have fewer stomatal conductance measurements than our other data types.

2.5 Isotopic analysis

We collected wood cores from 31 individual trees at a 1.5 m height using a 5.15 mm diameter increment borer (JIM-GEM, Forestry Suppliers Inc., Jackson, MS, USA). Since no definable tree rings were apparent, we created a fine powder for isotope analysis by drilling holes into dried cores using a dry ceramic drill bit (Dremel) along the outermost 5 cm of wood below the bark, which was chosen to represent the most recent carbon signal for ¹³C analyses. The fine powder (200 mesh, 0.2–5 mg) was then mixed, and a random sample was used to extract ¹³C∕¹²C ratios (to obtain δ¹³C values against the VPDB standard), which we estimated to be representative of at least the last 2–3 years, based on analogous literature growth rate values: O. xalapensis and A. acuminata range from 0.25 to 2.5 and from 0.6 to 0.9 cm yr⁻¹, respectively (Kappelle et al., 1996; Ortega-Pieck et al, 2011). These rates result in a 5 cm range of at least 2 and 5.5 years, though the high rates were determined for very young trees under very different conditions, and this is explicitly unknown in our study. Since we only sample the most recent years, no isotopic discrimination against atmospheric ¹³C due to preferential diffusion and carboxylation of ¹²C was conducted. Rather, we assume that δ¹³C values are representative of the relative amount of volcanic CO₂ vs. atmospheric CO₂ sequestered by the tree over the period of growth represented in the sample. δ¹³C values were determined by continuous-flow dual-isotope analysis using a CHNOS Elemental Analyzer and IsoPrime 100 mass spectrometer at the University of California Berkeley Center for Stable Isotope Biogeochemistry. External precision for C isotope determinations is ±0.10 ‰. Ten δ¹³C measurements did not have corresponding soil CO₂ flux measurements due to the flux measurements being unavailable for the final 2 days of sampling, and another five samples were from trees that showed signs of extreme stress, such as browning leaves or anomalously low fluorescence measurements. Since the purpose of our study was to explore the non-lethal effects of volcanic CO₂ on trees, during analysis we excluded all trees that were observed in the field as showing visible signs of stress or being not fully mature. After these exclusions, all remaining tree cores with co-located CO₂ flux measurements were from Turrialba.

2.6 Sulfur dioxide probability from satellite data

To assess the likelihood of trees having been significantly stressed in the past by volcanic sulfur dioxide (SO₂) from the central crater vents, we took two approaches. First, we were guided by in situ measurements taken in the same areas by Jenkins et al. (2012), who assessed the physiological interactions of SO₂ and CO₂ on vegetation on the upper slopes of Turrialba and demonstrated a rapid exponential decay of SO₂ away from the central vent. Second, for long-term exposure we derived the likelihood of exposure per unit area using satellite data sensitive to SO₂ (Fig. 2). The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), launched in December 1999 on NASA's Terra satellite, has bands sensitive to SO₂ emission in the thermal infrared (TIR), at ∼60 m × 60 m spatial resolution. We initially used ASTER Surface Radiance TIR data (AST_09T), using all ASTER observations of the target area over the entirety of the ASTER mission (October 2000 until writing began in late 2017). The TIR bands were corrected for downwelling sky irradiance and converted into units of W m⁻² µm⁻¹. For each observation, an absorption product is calculated by subtracting SO₂-insensitive from SO₂-sensitive bands:

$$\begin{array}{}\text{(1)}& S^t=\left(b_{\mathrm{10}}+b_{\mathrm{12}}\right)-\mathrm{2}\cdot b_{\mathrm{11}},\end{array}$$

where S is the SO₂ index, t is an index representing the time of acquisition, b₁₀ is the radiance at band 10 (8.125–8.475 µm), b₁₁ is the radiance at band 11 (8.475–8.825 µm), and b₁₂ is the radiance at band 12 (8.925–9.275 µm). This is similar to the method of Campion et al. (2010). The granules were then separated into day and night scenes, projected onto a common grid, and then thresholded to S > 0.1 W m⁻² µm⁻¹ and converted into a probability (Abrams et al., 2015). The output is a spatial dataset that describes the probability of an ASTER observation showing an absorption feature above a 0.1 W m⁻² µm⁻¹ threshold across the entirety of the ASTER observations for day or night separately. The number of scenes varies per target, but they tend to be between 200 and 800 observations in total, over the 17-year time period of satellite observations. However, certain permanent features, such as salt pans, show absorption features in band 11 and therefore have high ratios for the algorithm used. We therefore used a second method that seeks to map transient absorption features. For this method, we subtract the median from each S^t, yielding a median deviation stack. By plotting the maximum deviations across all observations, we then get a map of transient absorption features; in our case, these are mostly volcanic SO₂ plumes, which map out the cumulative position of different plume observations well. To speed up processing, some of the retrieval runs were binned in order to increase the signal-to-noise ratio, since the band difference can be rather noisy.

Figure 2 The influence of two potentially confounding gases on our study area (right-hand white shape) in Costa Rica is low to non-existent: anthropogenic CO₂ from San José (blue to red color scale), and volcanic SO₂ (purple color scale). White shapes are drawn around locations of the forested active volcanic edifices in Costa Rica. The red dashed line indicates the rough border of the San José urban area. Prevailing winds throughout the year consistently blow all anthropogenic CO₂ away from our study area and from all other white shapes.

2.7 Modeling the anthropogenic CO₂ influence from inventory data

We assessed the likelihood of anthropogenic CO₂, enhancements of air from San José, Costa Rica's capital and main industrial and population center, influencing our measurements. We used a widely applied FLEXible PARTicle dispersion model (Eckhardt et al., 2017; Stohl et al., 1998, 2005; Stohl and Thomson, 1999) in a forward mode (Stohl et al., 2005), FLEXPART, to simulate the downwind concentrations of CO₂ in the atmosphere (e.g., Belikov et al., 2016), due to inventory-derived fossil fuel (FF) emissions in our study area for the year 2015 (Fig. 2). The National Centers for Environmental Prediction (NCEP) – Climate Forecast System Reanalysis (CFSR) – 2.5^∘ horizontal resolution meteorology (Saha et al., 2010a, b) and 1 km Open-source Data Inventory for Anthropogenic CO₂ (ODIAC; Oda and Maksyutov, 2011) emissions for 2015 were used to drive the FlEXPART model. The CO₂ concentrations were generated at a 1 km spatial resolution within three vertical levels of the atmosphere (0–100, 100–300, and 300–500 m) that are possibly relevant to forest canopies in Costa Rica. However, to assess the magnitude of enhancements we only used CO₂ concentrations observed within the lowest modeled level of the atmosphere, from 0 to 100 m. Validation of the model with direct observations was not required because we were only interested in ensuring that anthropogenic CO₂ dispersed upslope from San José was not having a significant effect on our study area; we were not aiming to capture intra-canopy variability, typically present at tens to hundreds of parts per million variable, which is not relevant to the better mixed, distal single-digit-or-less parts per million signal from San José. The actual concentration of CO₂ and any biogenic influence in the modeled area was irrelevant because the spatial distribution of anthropogenic CO₂ was the only factor relevant for this test. 2015 was used as a representative year for simulating the seasonal cycle of CO₂ concentrations that would be present in any particular year.

3 Results

3.1 Volcanic CO₂ emissions through the soil

We measured CO₂ flux emitted through the soil at 66 points over 4 days (Fig. 1). The first eight points were on Irazú, and the rest were located near the Ariete fault on Turrialba. Mean soil CO₂ flux values over the entire sampling area varied from 3 to 37 g m⁻² d⁻¹, with an average of 11.6 g m⁻² d⁻¹ and a standard deviation of 6.6 g m⁻² d⁻¹. A 12-bin histogram of mean CO₂ flux shows a bimodal right-skewed distribution with a few distinct outliers (Fig. 3). Fluxes were generally larger on Irazú than on Turrialba. This result agrees with previous studies which showed that the north flank of Irazú has areas of extremely high degassing, whereas most of our sampling locations on Turrialba were in areas that had comparatively lower diffuse emissions (Epiard et al., 2017; Stine and Banks, 1991). We used a cumulative probability plot to identify different populations of CO₂ fluxes (Fig. 3; Cardellini et al., 2003; Sinclair, 1974).

Figure 3 Soil CO₂ flux into the sub-canopy air of forests on the Turrialba–Irazú volcanic complex is pervasively and significantly influenced by a deep volcanic gas source. At least four different overlapping populations of soil CO₂ flux were identified using a cumulative probability plot, where inflection points indicate population boundaries (Sinclair 1974). 69 % of sampling locations (45 total) are exposed to varying degrees of volcanically derived elevated CO₂. Populations are color-coded based on the same color scale as Fig. 1.

We created an inventory-based model of anthropogenic CO₂ emissions from the San José urban area, parts of which are less than 15 km from some of our sampling locations (Fig. 2). Our model shows that CO₂ emitted from San José is blown west to southwest by prevailing winds. Our study area is directly east of San José and as such is unaffected by anthropogenic CO₂ from San José, which is the only major urban area near Turrialba and Irazú. Since the trees sampled are spatially close to each other, they are exposed to the same regional background CO₂ variability. Additionally, we used ASTER data to map probabilities of SO₂ across Costa Rica as a possible confounding factor. The active craters of both Turrialba and Irazú emit measurable amounts of SO₂, which is reflected by the high SO₂ probabilities derived there (Fig. 2). Tropospheric SO₂ quickly converts to sulfate, a well-studied process intensified by the presence of volcanic mineral ash, plume turbulence, and a humid tropical environment (Oppenheimer et al., 1989; Eatough et al., 1994); furthermore, the bulk of the SO₂ emissions is carried aloft. Consequently, any remaining SO₂ causing acid damage effects on trees at Turrialba is limited to a narrow band of a few 100 m around the mostly quietly steaming central vent, which has been thoroughly ecologically evaluated for acid damage (Jenkins et al., 2012). D'Arcy (2018) has assessed this narrow, heavily SO₂-affected area immediately surrounding the central crater vent of Turrialba, which we avoided, and our sampling sites are mostly within their control zone that is not considered majorly affected by SO₂ but where diffuse CO₂ degassing dominates the excess gas phase (Epiard et al., 2017). Our study area is on the flanks of the volcano, where ASTER-derived SO₂ probability is minimal, and SO₂ influence not detectable on the ground (Jenkins et al., 2012; Campion et al., 2012). Most other volcanoes in Costa Rica emit little to no SO₂ on a decadal timescale, shown by the low or non-existent long-term SO₂probabilities over the other volcanoes in Costa Rica (white shapes in Fig. 2).

Figure 4 Bulk wood δ¹³C of trees on Costa Rica's Turrialba volcano shows strong correlations with increasing volcanic CO₂ flux for two species, O. xalapensis and B. nitida, indicating long-term photosynthetic incorporation of isotopically heavy volcanic CO₂. Stable carbon isotope ratio (δ¹³C) of wood cores are plotted against soil CO₂ flux measured immediately adjacent to the tree that the core sample was taken from. Background and volcanic influence labels apply to both axes – higher CO₂ flux and heavier (less negative) δ¹³C values are both characteristic of volcanic CO₂ emissions.

3.2 Tree core isotopes

Bulk wood δ¹³C measurements of all samples in this study, independent of exposure, ranged from −24.03 ‰ to −28.12 ‰, with most being clustered around −26 ‰ (Fig. 4). A five-bin histogram of all δ¹³C measurements shows a slightly right-skewed unimodal normal distribution, with an average of −26.37 ‰ and a standard deviation of 0.85 ‰. A. acuminata and O. xalapensis have nearly identical averages (−26.14 ‰ and −25.97 ‰, respectively), while B. nitida has a noticeably lighter average of −27.02 ‰. Diffuse excess CO₂ emissions throughout the investigation areas reflect a deep volcanic source which typically varies little in time (Epiard et al., 2017), but such diffuse emissions spatially follow geological subsurface structures (Giammanco et al., 1997). Their temporal variability therefore reflects long-term low-amplitude modulation of the volcanic heavy-δ¹³CO₂ signal, and their spatial distribution is mostly constant over tree lifetimes (Aiuppa et al., 2004; Peiffer et al., 2018; Werner et al., 2014), providing a constant long-term spatial gradient of CO₂ exposure to the forest canopy. Our data show that in areas where CO₂ flux is higher, the wood cores contained progressively higher amounts of ¹³C for two of the three species. Interestingly, our tree core δ¹³C showed no relationship with instantaneous stomatal conductance for any species, indicating that no stress threshold was exceeded during measurement across the sample set.

3.3 Plant function (fluorescence, chlorophyll, stomatal conductance)

Our measurements and literature data confirm that ecosystems growing in these locations are consistently exposed to excess volcanic CO₂, which may impact chlorophyll fluorescence, chlorophyll concentrations, and stomatal conductance of nearby trees. After excluding visibly damaged trees, leaf fluorescence, expressed as F_v∕F_m, was very high in most samples. F_v∕F_m ranged from 0.75 to 0.89, with most measurements clustering between 0.8 and 0.85 (Fig. 5). The fluorescence data have a left-skewed unimodal distribution. The leaf fluorescence (F_v∕F_m) values for A. acuminata had a strong positive correlation with soil CO₂ flux (r²=0.69, p<0.05), while the other two species showed no correlation. No confounding factors measured were correlated with F_v∕F_m for any species. In general, B. nitida had the highest F_v∕F_m values, and A. acuminata and O. xalapensis had similar values except for a few O. xalapensis outliers. Chlorophyll concentration measurements were highly variable, ranging from 260 to 922 µmol m⁻², with an average of 558 µmol m⁻² and a standard deviation of 162 µmol m⁻² (Fig. 6). Chlorophyll concentration had a complicated right-skewed bimodal distribution, likely due to the noticeably different averages for each species. A. acuminata and O. xalapensis both displayed weak correlations between chlorophyll concentration and soil CO₂ flux (r²=0.38 and r²=0.28, respectively), but their trend lines were found to be almost perpendicular (Fig. 6). As CO₂ flux increased, A. acuminata showed a slight increase in chlorophyll concentration, while O. xalapensis had significant decreases in chlorophyll concentration. B. nitida individuals growing on steeper slopes had significantly lower chlorophyll concentration measurements (r²=0.42, p<0.05) than those on gentler slopes, a trend not expressed by either of the other two species (r²=0.01 for both), demonstrating no significant influence of slope across the majority of samples. Stomatal conductance ranged from 83.5 to 361 mmol H₂O m⁻² s⁻¹, with an average of 214 mmol H₂O m⁻² s⁻¹ and a standard deviation of 73.5 mmol H₂O m⁻² s⁻¹. Distribution was bimodal, with peaks around 150 and 350 mmol H₂O m⁻² s⁻¹. A. acuminata had a moderate positive correlation (r²=0.51) with soil CO₂ flux, but it was not statistically significant due to a lack of data points (Fig. 7); however this is a result consistent with the observed higher chlorophyll concentration (Fig. 6). The other two species displayed no correlation with soil CO₂ flux. B. nitida had a moderate negative correlation (r²=0.61) with slope, similar to its correlation between chlorophyll concentration and slope.

Figure 5 Photosynthetic activity of some tree species in old-growth forests on the upper flanks of two active volcanoes in Costa Rica, Turrialba and Irazú, may show short-term response to volcanically elevated CO₂. Leaf fluorescence (F_v∕F_m) and soil CO₂ flux were strongly correlated for A. acuminata but not for other species.

Figure 6 Some tree species in old-growth forests on the upper flanks of two active volcanoes in Costa Rica, Turrialba and Irazú, may express their short-term response to volcanically elevated CO₂ by producing more chlorophyll. A species that showed strong short-term response (A. Acuminata; Fig. 5) also shows a positive correlation between chlorophyll concentration and mean soil CO₂ flux.

Figure 7 Leaf stomatal conductance of a tree species that strongly responds to volcanically elevated CO₂ (Figs. 5, 6) has positive correlations with volcanic CO₂ flux, consistent with increased gas exchange.

4 Discussion

4.1 Long-term plant uptake in volcanic CO₂

Turrialba and Irazú continuously emit CO₂ through their vegetated flanks, but prior to this study it was unknown if the trees growing there were utilizing this additional isotopically heavy volcanic CO₂. All tree cores with corresponding CO₂ flux measurements were from areas proximal to the Ariete fault on Turrialba, where atmospheric and volcanic δ¹³C have significantly different values (−9.2 ‰ and −3.4 ‰, respectively; Malowany et al., 2017). If the trees assimilate volcanic CO₂ through their stomata, then we would expect wood δ¹³C to trend towards heavier values as diffuse volcanic CO₂ flux increases. Studies at FACE sites have found that altering the isotopic composition of the air by artificially adding CO₂ with a different carbon isotope composition than the atmosphere leads to significant changes in the δ¹³C value of plant matter and tree rings growing there, leading us to expect similar effects from the naturally added volcanic CO₂ (Körner, 2005). It is worth noting that the FACE CO₂ (δ¹³C is −29.7 ‰), is significantly depleted in ¹³C compared the atmosphere, whereas volcanic CO₂ is enriched (δ¹³C is −3.4 ‰ at Turrialba) compared to atmosphere (Körner, 2005). After excluding damaged samples and stressed trees, δ¹³C was strongly correlated with soil CO₂ flux for both B. nitida and O. xalapensis (Fig. 4). A. acuminata did not have a statistically significant correlation between soil CO₂ flux and δ¹³C, likely because it had the fewest data points and a minimal range of CO₂ and δ¹³C values. The difference in regression slope between B. nitida and O. xalapensis (Fig. 4) may be due to physiological differences across traits or species, and/or due to differences in exposure owing to canopy height differences. Resolving this question would require a much larger multi-species sample size which could only be sufficiently obtained using remote sensing methods. The strong positive correlations between CO₂ flux and increasingly heavy δ¹³C values suggest that the trees have consistently photosynthesized with isotopically heavy excess volcanic CO₂ over the last few years and are therefore growing in eCO₂ conditions. Assuming that most of the variations in δ¹³C are caused by incorporation of heavy volcanic CO₂, we can calculate the average concentration of the mean volcanic excess CO₂ in the air the plants are exposed to with a mass balance equation (Eq. 2):

$$\begin{array}{}\text{(2)}& {\mathrm{C}}_{\mathrm{v}}={\displaystyle \frac{\left({\mathit{\delta }}_{\mathrm{b}}-{\mathit{\delta }}_t\right)}{({\mathit{\delta }}_{\mathrm{a}}-{\mathit{\delta }}_{\mathrm{v}})}}{\mathrm{C}}_{\mathrm{a}},\end{array}$$

where C_v is the mean volcanic excess component of the CO₂ concentration in air, C_a is the atmospheric “background” (i.e., non-volcanic) CO₂ concentration, δ_a is atmospheric δ¹³C, δ_b is the most negative δ¹³C measurement for the species being studied, δ_t is the δ¹³C value for the tree for which volcanic CO₂ exposure is being calculated, and δ_v is the δ¹³C value of the volcanic CO₂. Background wood δ¹³C is the value of the point for each species with the lowest CO₂ flux (Fig. 4), and the other wood δ¹³C measurement is any other point from the same species. Values for δ_v, δ_a, and C_a are taken from Malowany et al. (2017). For the tree core with the highest measured CO₂ flux for O. xalapensis, this equation yields a mean excess volcanic CO₂ concentration of 115 ppm, bringing the combined mean atmospheric (including volcanic) CO₂ concentration tree exposure to potentially around ∼520 ppm. For B. nitida this equation yields 133 ppm of mean excess volcanic CO₂ at the highest flux location, for a combined total mean of potentially ∼538 ppm CO₂. These numbers may be on the high side as the calculation assumes that carbon isotope discrimination remains constant for all trees within a given species, but they serve as an estimate of the approximate magnitude of the average amount of CO₂ that these trees are exposed to. A ¹⁴C tree ring study at Mammoth Mountain found an average yearly volcanic excess CO₂ exposure of 20–70 ppm over a 15-year period (Lewicki et al., 2014). Turrialba is significantly more active than Mammoth Mountain, so trees growing in high emission areas of Turrialba may be exposed to similar or higher amounts of CO₂ than the tree in the Mammoth Mountain study. Additional measurements of tree core δ¹³C and associated soil CO₂ fluxes would help corroborate our observations, which were based on a limited number of data points. Though tree ring ¹⁴C content in volcanically active areas has been linked to variations in volcanic CO₂ emissions, and comparing patterns of δ¹³C to ¹⁴C measurements for the same wood samples could provide additional confirmation of this finding (Evans et al., 2010; Lefevre et al., 2017; Lewicki et al., 2014), this additional dimension was outside the scope of this exploratory study. However, beyond such pattern confirmation, using ¹⁴C dating of trees exposed to naturally isotopically distinct excess CO₂ is, in fact, unfortunately not a reliable method for these environments due to the well-known δ¹⁴C deficiency in trees exposed to excess volcanic CO₂ which is isotopically “dead” with respect to ¹⁴C, creating spurious patterns that preclude dating by ¹⁴C (e.g., Lefevre et al., 2017; Lewicki et al., 2014).

Our data demonstrate that CO₂ fluxes through the soil may be a representative relative measure for eCO₂ exposure of overlying tree canopies. Forest canopy exposure to volcanic CO₂ will vary over time, as will volcanic eCO₂; once emitted through the soil into the sub-canopy atmosphere, the gas experiences highly variable thermal and wind disturbances which significantly affect dispersion of CO₂ on minute-to-minute, diurnal, and seasonal timescales (Staebler and Fitzjarrald, 2004; Thomas, 2011). These processes cause in-canopy measurements of CO₂ concentration to be highly variable, making instantaneous concentration measurements in a single field campaign not representative of long-term relative magnitudes of CO₂ exposure. Soil CO₂ fluxes are less tied to atmospheric conditions and are primarily externally modulated by rainfall which increases soil moisture and therefore lowers the soil's gas permeability (Camarda et al., 2006; Viveiros et al., 2009). These fluxes can also be affected by variations in barometric pressure, but both of these factors are easily measurable and therefore can be factored in when conducting field work (Viveiros et al., 2009). Assuming the avoidance of significant rainfall and pressure spikes during sampling (measurements were conducted in the dry season and no heavy rains or significant meteorological variations in pressure occurred during field work), measuring the input of CO₂ into the sub-canopy atmosphere as soil CO₂ fluxes is therefore expected to better represent long-term input and exposure of tree canopies to eCO₂ than direct instantaneous measurements of sub-canopy CO₂ concentration. Previous studies at Turrialba have shown that local volcanic CO₂ flux is relatively constant on monthly to yearly timescales (de Moor et al., 2016). Therefore, current soil CO₂ fluxes should give relatively accurate estimates of CO₂ exposure over time. This paper corroborates that expectation by demonstrating strong spatial correlations between volcanically enhanced soil CO₂ emissions with co-located stable carbon isotope signals of these emissions documented in the trees' xylem.

A study at the previously mentioned Mammoth Mountain tree kill area examined the connection between δ¹³C and volcanic CO₂ fluxes but focused on the difference between trees killed by extreme CO₂ conditions and those that were still alive (Biondi and Fessenden, 1999). They concluded that the changes in δ¹³C that they observed were due to extreme concentrations of CO₂ (soil CO₂ concentrations of up to 100 %) impairing the functioning of root systems, leading to closure of stomata and water stress (Biondi and Fessenden, 1999). CO₂ does not inherently harm trees, but the extreme CO₂ concentrations (up to 100 % soil CO₂) at the Mammoth Mountain area caused major soil acidification, which led to the tree kill (McGee and Gerlach, 1998). We have evidence that those acidification processes are not affecting our δ¹³C measurements and that variations in our δ¹³C measurements are more likely to be caused by the direct photosynthetic incorporation of isotopically heavy volcanic CO₂. Our δ¹³C measurements have no statistically significant correlation with stomatal conductance, which suggests that our heavier δ¹³C measurements are not linked to stomatal closure. None of the trees included in the analysis displayed obvious signs of stress, from water or other factors, as indicated by their high fluorescence and chlorophyll concentration values and the lack of visible indicators of stress; specifically, our values of F_v∕F_m ∼0.8 indicate that the PSII was operating efficiently in most of the trees we measured (Baker and Oxborough, 2004). The Mammoth Mountain tree kill areas have higher CO₂ fluxes by several orders of magnitude (well over 10 000 g m⁻² d⁻¹) than the areas we sampled (up to 38 g m⁻² d⁻¹), making it much more likely that stress from soil acidification causes stomatal closure and affects wood δ¹³C measurements at Mammoth Mountain (Biondi and Fessenden, 1999; McGee and Gerlach, 1998; Werner et al., 2014). In contrast, most of the diffuse degassing at Turrialba does not lead to soil acidification or pore space saturation, as is evident in our own and others' field data (e.g., Epiard et al., 2017). Thus, changes in our δ¹³C values are best explained by direct photosynthetic incorporation of isotopically heavy volcanic CO₂. To the best of our knowledge, this is the first time that a direct correlation between volcanic soil CO₂ flux and wood δ¹³C has been documented. Future studies should explore this correlation further, as our findings are based on a limited sample size.

4.2 Short-term species response to eCO₂

Short-term plant functional responses at the leaf level to elevated CO₂ were highly species-dependent. B. nitida had no statistically significant functional responses to soil CO₂ flux, and O. xalapensis only had a weak negative correlation between soil CO₂ flux and chlorophyll concentration (Fig. 6). A. acuminata, a nitrogen-fixing species, was the only species with a consistent and positive functional response to elevated CO₂, displaying a strong positive correlation with fluorescence and a weak positive correlation with chlorophyll concentration and stomatal conductance (Figs. 5–7). Previous studies which linked changes in the NDVI to pre-eruptive volcanic activity on the flanks of Mt. Etna and Mt. Nyiragongo support our observation of a correlation between plant function and volcanic CO₂ flux (Houlié et al., 2006; Seiler et al., 2017). This link raises the question of why only one of three species displayed strong functional responses to volcanic CO₂. The lack of response in B. nitida and O. xalapensis could be due to nitrogen limitation, a factor that would not affect A. acuminata due to its nitrogen-fixing capability. Previous studies have found that nitrogen availability strongly controls plant responses to both naturally and artificially elevated CO₂ concentrations in a variety of ecosystems, including grasslands and temperate forests (Garten et al., 2011; Hebeisen et al., 1997; Lüscher et al., 2000; Norby et al., 2010; Tognetti et al., 2000). Nitrogen limitation has been posited to be an important factor in tropical montane cloud forests and may be contributing to the lack of responses in B. nitida and O. xalapensis (Tanner et al., 1998). Due to the exploratory nature of our study, we do not have a large enough dataset to conclude that the nitrogen-fixing capability of species like A. acuminata is the cause for its positive response to volcanically elevated CO₂ concentrations, as has been speculated before (Schwandner et al., 2004), but it is a possible correlation that deserves further investigation.

4.3 Time constraints

To support these results, we further assessed the possibility of effects of time constraints on growth rates and isotopic signals, despite the compelling spatial variability in the independent variable (naturally isotopically labeled excess volcanic CO₂) in our study (Helle and Schleser, 2004; Verheyden et al., 2004). As tropical trees typically lack tree rings, it is difficult to directly constrain the precise time period that the data represent. However, since we sampled from the outside in, all the samples appear to at least have the most recent growth period in common. To assess how far back in time our samples could likely represent, we compared our sampled core depths to reported growth rates for the same species in similar environments. Reported growth rates for two of our species, O. xalapensis and A. acuminata, range from 0.25 to 2.5 and 0.6 to 0.9 cm yr⁻¹, respectively (Kappelle et al., 1996; Ortega-Pieck et al., 2011). Given that our samples are bulk measurements of the outer 5 cm of wood, each sample would represent between 2 and 5.5 years, although the conditions that these growth rates were measured in were different than in our study. Clear time constraints would be necessary for higher resolution analysis, but this need is somewhat mitigated by the continuous, long-term, and, over multiple decades, mostly invariant nature of diffuse volcanic CO₂ emissions, which is completely independent of any non-volcanic environmental influences on growth rates. By providing an upper and lower bound in the expected growth span represented in our samples, we believe that these samples represent similar time frames during the continuous exposure to excess volcanic CO₂ over the lifetimes of the trees sampled. Due to the continuous nature of the volcanic CO₂ enhancement, we are not investigating and analyzing transient events, and our results instead represent spatial variability in excess CO₂ availability averaged over similar time periods.

Although we do not believe that our samples represent a long enough time period for long-term variations in δ¹³C (Suess effect) to be relevant, if it does affect our samples, it would be beneficial for detection of volcanic CO₂, as the Suess effect is gradually increasing the gap between atmospheric and volcanic δ¹³C. Since our δ¹³C values likely represent several years of growth, small-scale temporal variations in excess volcanic CO₂ release are unlikely to significantly impact the results. Larger trees tend to grow slower than smaller trees, so the outer 5 cm of wood should represent a longer time period in larger trees. Thus, if temporal variations had a significant effect on our δ¹³C measurements, we would expect this to be represented by some correlation between DBH and δ¹³C, which is not present for any species studied. Three of the five B. nitida individuals measured were very large (150–190 cm DBH), whereas the other two are much smaller (11.5 and 15.3 cm DBH). Although the age and growth rates of these two groups of trees likely vary significantly, we found no correlation between DBH and δ¹³C, though we did find a strong correlation between the completely independent diffuse excess (volcanic) CO₂ flux and wood δ¹³C. Furthermore, the relationships presented are on a per species basis to avoid complications resulting from different growth rates across species. This is important because δ¹³C values provide an integral value of assimilated carbon by the entire tree (not just individual leaves). The depth of tree core sample was identical for each species (the outermost part of the trunk), and we can safely assume that the volcanic CO₂ exposure has been consistent over the time period under investigation.

Because individual time variability of growth rates can possibly affect these signals as well, future studies that attempt to study tree ring isotopes in this context at higher resolutions will likely require stricter and more detailed time constraints and cell-level stress analysis to average out the effects of long-term variations in δ¹³C (Suess effect), seasonal cycles, potential short-term transient stress-induced growth rate variations, effects of water use efficiency (WUE), and potential short-term variations in CO₂ flux, all of which may result in time-averaged isotopic shifts over different growth periods (Helle and Schleser, 2004; Verheyden et al., 2004). We include these notes as guidance in Sect. 4.4: Lessons learned for future studies. Despite the additional difficulty of conducting higher time resolution analysis, this type of study holds great potential for attempting to reconstruct volcanic CO₂ histories and to study its potential fertilization effect, due to the completely independent nature of the volcanic excess CO₂ supply to the sub-canopy air.

4.4 Lessons learned for future studies

This exploratory study reveals significant new potential for future studies to utilize the volcanically enhanced CO₂ emissions approach to study tropical ecosystem responses to eCO₂ – one of the largest uncertainties in climate projections. Costa Rica's volcanoes are host to large areas of relatively undisturbed rainforest, making them ideal study areas for examining responses of ecosystems to eCO₂. However, there are several challenges that future studies should take into consideration if attempting to expand upon this preliminary study. Given the enormous tropical species diversity and the need to control for confounding factors, large datasets will be needed to answer these questions conclusively. One open question, for example, is how WUE in upper and lower canopy leaves of the same and different individuals within a species may affect isotopic sequestration of CO₂. Since the excess volcanic CO₂ is naturally isotopically labeled, this could be assessed by a much more detailed (by individual tree leaf, branch, and xylem) core study coupled with long-term measurements of evapotranspiration, heat stress, and stomatal conductance, the last of which in our study showed no significant correlation with the δ¹³C signal in the wood cores across spatial gradients. Field data can be difficult to acquire in these rugged and challenging environments. A remote sensing approach using airborne measurements, validated by targeted representative ground campaigns, could provide sufficiently large datasets to represent species diversity and conditions appropriately. Many of the data types that would be useful for this type of study can be acquired from airborne platforms, and remote sensing instruments can quickly produce the massive datasets required to provide more comprehensive answers to these questions. A recent meta-analysis showed that studies at natural CO₂-producing springs and FACE experiments have found similar results in a variety of plant traits, which significantly strengthens the case that volcanoes are a potentially extremely valuable resource for determining plant responses to elevated CO₂ concentrations (Saban et al., 2019). While the spring studies have yielded valuable results, volcanoes could offer several advantages over springs for future studies. Active volcanoes are significantly larger systems than non-volcanic springs and often feature several CO₂-producing springs and also several dry gas seeps; this offers more data, more control points for comparison, greater species diversity, and greater potential for comprehensive measurements of a statistically meaningful dataset from remote sensing platforms. Due to their volcanic hazard potential, volcanoes are also more likely to already have long-term volcanological monitoring programs for CO₂ fluxes and ecological disturbances, which may be utilized to analyze the long-term effects of enhanced levels of CO₂ emissions on these volcanically active tropical ecosystems.

Our results also offer significant new tools for the volcanology, where reconstructing past volcano behavior through eruption histories is hampered by severe preservation gaps in the stratigraphic record. A strong link between δ¹³C and volcanic CO₂ could be a game changer by establishing long-term histories of volcanic CO₂ emission variations. These proxy signals could be traced back in time using living and preserved dead trees in order to fill gaps in the historical and monitoring records – a boon for volcano researchers and observatories to improve eruption prediction capabilities (Newhall et al., 2017; Pyle, 2017; Sparks et al., 2012). While variations in tree ring ¹⁴C content have been shown to correlate well with variations in volcanic CO₂ flux (Evans et al., 2010; Lefevre et al., 2017; Lewicki and Hilley, 2014), ¹⁴C is relatively expensive to measure, limiting the spatial and temporal coverage of data that can be acquired. ¹³C is an inexpensive alternative to ¹⁴C and can be measured at more laboratories, allowing for substantially more data to be acquired. Some previously mentioned studies (Lefevre et al., 2017; Pasquier-Cardin et al., 1999) have found correlations between ¹³C and ¹⁴C in plants that have incorporated volcanic CO₂, strengthening the potential for using ¹³C in this type of study. Further development of the ¹³C approach to tracking volcanic CO₂ emissions would prove beneficial to future studies attempting to use plants to study large areas and timescales of volcanic degassing. Independent validation, and calibration by wood core dendrochronology via ¹⁴C, tree rings, or chemical event tracers like sulfur isotopes, could significantly advance the concept of using wood carbon as an archive of past degassing activity. Crucially, these tree ring archives could provide temporal records of degassing at dangerous volcanoes which have previously been poorly monitored or not monitored at all, significantly improving the accuracy of hazard assessments. Furthermore, knowledge of the short-term real-time response of leaves to diffusely emitted eCO₂, which is more likely to represent deeper processes inside volcanoes than crater-area degassing (Camarda et al., 2012), may permit the use of trees as sensors of transient changes in volcanic degassing indicative of volcanic reactivation and deep magma movement possibly leading up to eruptions (Camarda et al., 2012; Houlié et al., 2006; Pieri et al., 2016; Schwandner et al., 2017; Seiler et al., 2017; Shinohara et al., 2008; Werner et al., 2013). To the best of our knowledge, we are the first to propose utilizing the combination of short-term leaf functional responses to volcanic CO₂ with long-term changes in δ¹³C values of wood for assessment of past and present volcanic activity in a single study.

5 Conclusions

Multiple areas of dense tropical forest on two Costa Rican active volcanoes are consistently and continuously exposed to volcanically elevated levels of atmospheric CO₂, diffusively cold emitted through soils into overlying forests. These isotopically heavy volcanic CO₂ emissions, which are mostly invariant, not accompanied by acidic gases, and independent of processes affecting growth rates, are well correlated with increases in heavy carbon signatures in wood cores from two species of tropical trees, possibly suggesting long-term incorporation of enhanced levels of volcanically emitted CO₂ into biomass. Each tree studied was co-located with a soil CO₂ flux measurement, and their soil CO₂ flux signals vary spatially around a continuous long-term local natural excess volcanic CO₂ source, which creates a local CO₂ gradient within which all the sampled trees are found. The excess volcanic CO₂ through local fault-bound gas seeps provides continuous exposure to all sampled trees over timescales much greater than the lifetimes of individual trees. Based on our limited exploratory measurements, confounding factors that are known to influence δ¹³C values in wood appear not to have significantly affected our measurements, indicating that the heavier wood δ¹³C values could be caused by photosynthetic incorporation of volcanic excess CO₂. One of the three species studied (A. acuminata) has consistent positive correlations between instantaneous plant function measurements and diffuse CO₂ flux measurements, indicating that short-term variations in elevated CO₂ emissions may measurably affect trees growing in areas of diffuse volcanic gas emissions. These observations reveal significant potential for future studies to use these areas of naturally elevated CO₂ to study ecosystem responses to elevated CO₂ and to use trees as sensors of the changing degassing behavior of volcanic flanks, which is indicative of deep magmatic processes.