Plant responses to volcanically elevated CO2 in two Costa Rican forests

Abstract

We explore the use of active volcanoes to determine the short- and long-term effects of elevated ${\mathrm{CO}}_{\mathrm{2}}$ on tropical trees. Active volcanoes continuously but variably emit ${\mathrm{CO}}_{\mathrm{2}}$ through diffuse emissions on their flanks, exposing the overlying ecosystems to elevated levels of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$. We found tight correlations ($r^{\mathrm{2}}=\mathrm{0.86}$ and $r^{\mathrm{2}}=\mathrm{0.74}$) between wood stable carbon isotopic composition and co-located volcanogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions for two of three investigated species (Oreopanax xalapensis andBuddleja 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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${\mathrm{CO}}_{\mathrm{2}}$ in the tropics and elsewhere. Results also point the way toward a possible future utilization of ecosystems exposed to volcanically elevated${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ tropical ecosystems can absorb; this is known as the “${\mathrm{CO}}_{\mathrm{2}}$ fertilization effect” (Gregory et al., 2009; Kauwe et al., 2016; Keeling, 1973; Schimel et al., 2015).

Free Air ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$-emitting natural springs to study plant responses to elevated ${\mathrm{CO}}_{\mathrm{2}}$, but these have been limited in scope due to the small spatial areas around springs that experience elevated ${\mathrm{CO}}_{\mathrm{2}}$ (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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ emissions made control points more available (Sharma and Williams, 2009; Tercek et al., 2008). These studies reported changes in rubisco, an enzyme central to ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ on ecosystems (Cawse-Nicholson et al., 2018; Vodnik et al., 2018). The broad flanks of active volcanoes experience diffuse emissions of excess ${\mathrm{CO}}_{\mathrm{2}}$ because the underlying active magma bodies continuously release gas, dominated by ${\mathrm{CO}}_{\mathrm{2}}$ 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${}^{\mathrm{2}}$, 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., ${\mathrm{SO}}_{\mathrm{2}}$, HF, HCl, and ${\mathrm{H}}_{\mathrm{2}}$S) have been scrubbed out in the deep subsurface, leading to a diffusely emanated gas mix of predominantly ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ (${\mathrm{eCO}}_{\mathrm{2}}$), although the specific effects of this ${\mathrm{eCO}}_{\mathrm{2}}$ on the trees are not well understood. Volcanic ${\mathrm{CO}}_{\mathrm{2}}$ has no ${}^{\mathrm{14}}\mathrm{C}$ and a ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ signature typically ranging from around $-$7 ‰ to $-$1 ‰, which is distinct from typical vegetation and noticeably enriched in ${}^{\mathrm{13}}\mathrm{C}$ compared to typical atmospheric values (Mason et al., 2017). If trees incorporate volcanic ${\mathrm{CO}}_{\mathrm{2}}$, then the stable carbon isotopic composition of wood may document the long-term, possibly variable influence of volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ on ecosystems, especially in understudied regions like the tropics. Several studies have found correlations between variations in volcanic ${\mathrm{CO}}_{\mathrm{2}}$ flux and plant ${}^{\mathrm{14}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ flux to decreases in ${}^{\mathrm{14}}\mathrm{C}$ content in tree rings in one or two trees, demonstrating the method's utility for uncovering yearly scale variations in volcanic ${\mathrm{CO}}_{\mathrm{2}}$ fluxes (Cook et al., 2001; Evans et al., 2010; Lewicki et al., 2014). The Naples study instead focused on using ${}^{\mathrm{14}}\mathrm{C}$ in grasses as short-term (2 to 6 months) monitors of volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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^{\mathrm{2}}\mathit{>}\mathrm{0.85}$) correlations between depletions in ${}^{\mathrm{14}}\mathrm{C}$ and enrichments in ${}^{\mathrm{13}}\mathrm{C}$ from volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${}^{\mathrm{13}}\mathrm{C}$ and ${}^{\mathrm{14}}\mathrm{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 ${\mathrm{eCO}}_{\mathrm{2}}$ might be revealed by plant functional measurements at the leaf scale, where the additional ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ flux across north flank of Irazú (a) and south flank of Turrialba (b). Colors of dots correspond to flux populations (see Fig. 3).

[Figure omitted. See PDF]

2 Methods

2.1 Investigated locations and sampling strategy

Irazú and Turrialba are two active volcanoes located $\sim \mathrm{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${}^{\mathrm{2}}$. 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 ${\mathrm{CO}}_{\mathrm{2}}$ emissions in areas proximal to the crater were calculated at $\mathrm{113}\pm \mathrm{46}$ tons d${}^{-\mathrm{1}}$ (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 ${\mathrm{CO}}_{\mathrm{2}}$ emissions on Turrialba (Epiard et al., 2017; Rizzo et al., 2016). Atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ has an average ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ value of $-$9.2 ‰ at Turrialba, and the volcanic ${\mathrm{CO}}_{\mathrm{2}}$ released at the Ariete fault has significantly heavier ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ represent the vast majority of gas discharge from Irazú, as the main crater releases 3.8 t d${}^{-\mathrm{1}}$ of ${\mathrm{CO}}_{\mathrm{2}}$ and a small area on the north flank alone releases 15 t d${}^{-\mathrm{1}}$ (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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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” ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{SO}}_{\mathrm{2}}$ concentrations, and local anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ from farms or that were likely to have heavily acidified soil. Excessively high soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ concentrations rather than direct consequences of elevated ${\mathrm{CO}}_{\mathrm{2}}$ on plants, we avoided areas with ${\mathrm{CO}}_{\mathrm{2}}$ fluxes high enough to possibly cause noticeable ${\mathrm{CO}}_{\mathrm{2}}$-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

${\mathrm{CO}}_{\mathrm{2}}$ concentrations and soil diffuse flux measurements

Soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ flux value to compare to the ${}^{\mathrm{13}}\mathrm{C}$ measurement of the corresponding tree sample. This technique is intended to provide a simple relative way to compare the ${\mathrm{CO}}_{\mathrm{2}}$ exposure of different trees, as a tree with high ${\mathrm{CO}}_{\mathrm{2}}$ flux near its base should experience consistently higher ${\mathrm{CO}}_{\mathrm{2}}$ concentrations than a tree with lower ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ fluxes, not concentrations, because the diffuse emissions of excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ sensor (LI-COR Inc., Lincoln, NE, USA) to measure soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$. 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_{\mathrm{o}}$) under conditions where we assume that the photosystem II (PSII) was fully reduced. Immediately following the $F_{\mathrm{o}}$ measurement, a 6000 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ saturation pulse was delivered from an array of red LEDs at 660 nm to record maximal fluorescence emission ($F_{\mathrm{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_{\mathrm{v}}/F_{\mathrm{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_{\mathrm{v}}/F_{\mathrm{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 ($\mathrm{\mu }$mol m${}^{-\mathrm{2}}$) 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_{\mathrm{s}}$ (mmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$), 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 ${}^{\mathrm{13}}\mathrm{C}$ analyses. The fine powder (200 mesh, 0.2–5 mg) was then mixed, and a random sample was used to extract ${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratios (to obtain ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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${}^{-\mathrm{1}}$, 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 ${}^{\mathrm{13}}\mathrm{C}$ due to preferential diffusion and carboxylation of ${}^{\mathrm{12}}\mathrm{C}$ was conducted. Rather, we assume that ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values are representative of the relative amount of volcanic ${\mathrm{CO}}_{\mathrm{2}}$ vs. atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ sequestered by the tree over the period of growth represented in the sample. ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 $\pm \mathrm{0.10}$ ‰. Ten ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurements did not have corresponding soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 (${\mathrm{SO}}_{\mathrm{2}}$) 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 ${\mathrm{SO}}_{\mathrm{2}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ on vegetation on the upper slopes of Turrialba and demonstrated a rapid exponential decay of ${\mathrm{SO}}_{\mathrm{2}}$ away from the central vent. Second, for long-term exposure we derived the likelihood of exposure per unit area using satellite data sensitive to ${\mathrm{SO}}_{\mathrm{2}}$ (Fig. 2). The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), launched in December 1999 on NASA's Terra satellite, has bands sensitive to ${\mathrm{SO}}_{\mathrm{2}}$ emission in the thermal infrared (TIR), at $\sim \mathrm{60}$ m $\times$ 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${}^{-\mathrm{2}}$ $\mathrm{\mu }$m${}^{-\mathrm{1}}$. For each observation, an absorption product is calculated by subtracting ${\mathrm{SO}}_{\mathrm{2}}$-insensitive from ${\mathrm{SO}}_{\mathrm{2}}$-sensitive bands:

1$$S^t=\left(b_{\mathrm{10}}+b_{\mathrm{12}}\right)-\mathrm{2}\cdot b_{\mathrm{11}},$$ where $S$ is the ${\mathrm{SO}}_{\mathrm{2}}$ index, $t$ is an index representing the time of acquisition, $b_{\mathrm{10}}$ is the radiance at band 10 (8.125–8.475 $\mathrm{\mu }$m), $b_{\mathrm{11}}$ is the radiance at band 11 (8.475–8.825 $\mathrm{\mu }$m), and $b_{\mathrm{12}}$ is the radiance at band 12 (8.925–9.275 $\mathrm{\mu }$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>\mathrm{0.1}$ W m${}^{-\mathrm{2}}$ $\mathrm{\mu }$m${}^{-\mathrm{1}}$ 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${}^{-\mathrm{2}}$ $\mathrm{\mu }$m${}^{-\mathrm{1}}$ 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 ${\mathrm{SO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ from San José (blue to red color scale), and volcanic ${\mathrm{SO}}_{\mathrm{2}}$ (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 ${\mathrm{CO}}_{\mathrm{2}}$ away from our study area and from all other white shapes.

[Figure omitted. See PDF]

2.7

Modeling the anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ influence from inventory data

We assessed the likelihood of anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$, 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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${}^{\circ }$ horizontal resolution meteorology (Saha et al., 2010a, b) and 1 km Open-source Data Inventory for Anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ (ODIAC; Oda and Maksyutov, 2011) emissions for 2015 were used to drive the FlEXPART model. The ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ and any biogenic influence in the modeled area was irrelevant because the spatial distribution of anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ was the only factor relevant for this test. 2015 was used as a representative year for simulating the seasonal cycle of ${\mathrm{CO}}_{\mathrm{2}}$ concentrations that would be present in any particular year.

3 Results

3.1

Volcanic ${\mathrm{CO}}_{\mathrm{2}}$ emissions through the soil

We measured ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ flux values over the entire sampling area varied from 3 to 37 g m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$, with an average of 11.6 g m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ and a standard deviation of 6.6 g m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$. A 12-bin histogram of mean ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ fluxes (Fig. 3; Cardellini et al., 2003; Sinclair, 1974).

Figure 3

Soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$. Populations are color-coded based on the same color scale as Fig. 1.

[Figure omitted. See PDF]

We created an inventory-based model of anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ variability. Additionally, we used ASTER data to map probabilities of ${\mathrm{SO}}_{\mathrm{2}}$ across Costa Rica as a possible confounding factor. The active craters of both Turrialba and Irazú emit measurable amounts of ${\mathrm{SO}}_{\mathrm{2}}$, which is reflected by the high ${\mathrm{SO}}_{\mathrm{2}}$ probabilities derived there (Fig. 2). Tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ 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 ${\mathrm{SO}}_{\mathrm{2}}$ emissions is carried aloft. Consequently, any remaining ${\mathrm{SO}}_{\mathrm{2}}$ 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 ${\mathrm{SO}}_{\mathrm{2}}$-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 ${\mathrm{SO}}_{\mathrm{2}}$ but where diffuse ${\mathrm{CO}}_{\mathrm{2}}$ degassing dominates the excess gas phase (Epiard et al., 2017). Our study area is on the flanks of the volcano, where ASTER-derived ${\mathrm{SO}}_{\mathrm{2}}$ probability is minimal, and ${\mathrm{SO}}_{\mathrm{2}}$ influence not detectable on the ground (Jenkins et al., 2012; Campion et al., 2012). Most other volcanoes in Costa Rica emit little to no ${\mathrm{SO}}_{\mathrm{2}}$ on a decadal timescale, shown by the low or non-existent long-term SO${}_{\mathrm{2}}$probabilities over the other volcanoes in Costa Rica (white shapes in Fig. 2).

Figure 4

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

[Figure omitted. See PDF]

3.2 Tree core isotopes

Bulk wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ 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-${\mathit{\delta }}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ exposure to the forest canopy. Our data show that in areas where ${\mathrm{CO}}_{\mathrm{2}}$ flux is higher, the wood cores contained progressively higher amounts of ${}^{\mathrm{13}}\mathrm{C}$ for two of the three species. Interestingly, our tree core ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$, which may impact chlorophyll fluorescence, chlorophyll concentrations, and stomatal conductance of nearby trees. After excluding visibly damaged trees, leaf fluorescence, expressed as $F_{\mathrm{v}}/F_{\mathrm{m}}$, was very high in most samples. $F_{\mathrm{v}}/F_{\mathrm{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_{\mathrm{v}}/F_{\mathrm{m}}$) values for A. acuminata had a strong positive correlation with soil ${\mathrm{CO}}_{\mathrm{2}}$ flux ($r^{\mathrm{2}}=\mathrm{0.69}$, $p\mathit{<}\mathrm{0.05}$), while the other two species showed no correlation. No confounding factors measured were correlated with $F_{\mathrm{v}}/F_{\mathrm{m}}$ for any species. In general, B. nitida had the highest $F_{\mathrm{v}}/F_{\mathrm{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 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$, with an average of 558 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ and a standard deviation of 162 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ (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 ${\mathrm{CO}}_{\mathrm{2}}$ flux ($r^{\mathrm{2}}=\mathrm{0.38}$ and $r^{\mathrm{2}}=\mathrm{0.28}$, respectively), but their trend lines were found to be almost perpendicular (Fig. 6). As ${\mathrm{CO}}_{\mathrm{2}}$ 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^{\mathrm{2}}=\mathrm{0.42}$, $p\mathit{<}\mathrm{0.05}$) than those on gentler slopes, a trend not expressed by either of the other two species ($r^{\mathrm{2}}=\mathrm{0.01}$ for both), demonstrating no significant influence of slope across the majority of samples. Stomatal conductance ranged from 83.5 to 361 mmol ${\mathrm{H}}_{\mathrm{2}}$O m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, with an average of 214 mmol ${\mathrm{H}}_{\mathrm{2}}$O m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ and a standard deviation of 73.5 mmol ${\mathrm{H}}_{\mathrm{2}}$O m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. Distribution was bimodal, with peaks around 150 and 350 mmol ${\mathrm{H}}_{\mathrm{2}}$O m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. A. acuminata had a moderate positive correlation ($r^{\mathrm{2}}=\mathrm{0.51}$) with soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ flux. B. nitida had a moderate negative correlation ($r^{\mathrm{2}}=\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$. Leaf fluorescence ($F_{\mathrm{v}}/F_{\mathrm{m}}$) and soil ${\mathrm{CO}}_{\mathrm{2}}$ flux were strongly correlated for A. acuminata but not for other species.

[Figure omitted. See PDF]

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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ flux.

[Figure omitted. See PDF]

Figure 7

Leaf stomatal conductance of a tree species that strongly responds to volcanically elevated ${\mathrm{CO}}_{\mathrm{2}}$ (Figs. 5, 6) has positive correlations with volcanic ${\mathrm{CO}}_{\mathrm{2}}$ flux, consistent with increased gas exchange.

[Figure omitted. See PDF]

4 Discussion

4.1

Long-term plant uptake in volcanic ${\mathrm{CO}}_{\mathrm{2}}$

Turrialba and Irazú continuously emit ${\mathrm{CO}}_{\mathrm{2}}$ through their vegetated flanks, but prior to this study it was unknown if the trees growing there were utilizing this additional isotopically heavy volcanic ${\mathrm{CO}}_{\mathrm{2}}$. All tree cores with corresponding ${\mathrm{CO}}_{\mathrm{2}}$ flux measurements were from areas proximal to the Ariete fault on Turrialba, where atmospheric and volcanic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ have significantly different values ($-$9.2 ‰ and $-$3.4 ‰, respectively; Malowany et al., 2017). If the trees assimilate volcanic ${\mathrm{CO}}_{\mathrm{2}}$ through their stomata, then we would expect wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ to trend towards heavier values as diffuse volcanic ${\mathrm{CO}}_{\mathrm{2}}$ flux increases. Studies at FACE sites have found that altering the isotopic composition of the air by artificially adding ${\mathrm{CO}}_{\mathrm{2}}$ with a different carbon isotope composition than the atmosphere leads to significant changes in the ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ value of plant matter and tree rings growing there, leading us to expect similar effects from the naturally added volcanic ${\mathrm{CO}}_{\mathrm{2}}$ (Körner, 2005). It is worth noting that the FACE ${\mathrm{CO}}_{\mathrm{2}}$ (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ is $-$29.7 ‰), is significantly depleted in ${}^{\mathrm{13}}\mathrm{C}$ compared the atmosphere, whereas volcanic ${\mathrm{CO}}_{\mathrm{2}}$ is enriched (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ is $-$3.4 ‰ at Turrialba) compared to atmosphere (Körner, 2005). After excluding damaged samples and stressed trees, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ was strongly correlated with soil ${\mathrm{CO}}_{\mathrm{2}}$ flux for both B. nitida and O. xalapensis (Fig. 4). A. acuminata did not have a statistically significant correlation between soil ${\mathrm{CO}}_{\mathrm{2}}$ flux and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, likely because it had the fewest data points and a minimal range of ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ flux and increasingly heavy ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values suggest that the trees have consistently photosynthesized with isotopically heavy excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ over the last few years and are therefore growing in ${\mathrm{eCO}}_{\mathrm{2}}$ conditions. Assuming that most of the variations in ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ are caused by incorporation of heavy volcanic ${\mathrm{CO}}_{\mathrm{2}}$, we can calculate the average concentration of the mean volcanic excess ${\mathrm{CO}}_{\mathrm{2}}$ in the air the plants are exposed to with a mass balance equation (Eq. 2): 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}},$$ where $C_{\mathrm{v}}$ is the mean volcanic excess component of the ${\mathrm{CO}}_{\mathrm{2}}$ concentration in air, C${}_{\mathrm{a}}$ is the atmospheric “background” (i.e., non-volcanic) ${\mathrm{CO}}_{\mathrm{2}}$ concentration, ${\mathit{\delta }}_{\mathrm{a}}$ is atmospheric ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, ${\mathit{\delta }}_{\mathrm{b}}$ is the most negative ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurement for the species being studied, ${\mathit{\delta }}_{\mathrm{t}}$ is the ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ value for the tree for which volcanic ${\mathrm{CO}}_{\mathrm{2}}$ exposure is being calculated, and ${\mathit{\delta }}_{\mathrm{v}}$ is the ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ value of the volcanic ${\mathrm{CO}}_{\mathrm{2}}$. Background wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ is the value of the point for each species with the lowest ${\mathrm{CO}}_{\mathrm{2}}$ flux (Fig. 4), and the other wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurement is any other point from the same species. Values for ${\mathit{\delta }}_{\mathrm{v}}$, ${\mathit{\delta }}_{\mathrm{a}}$, and C${}_{\mathrm{a}}$ are taken from Malowany et al. (2017). For the tree core with the highest measured ${\mathrm{CO}}_{\mathrm{2}}$ flux for O. xalapensis, this equation yields a mean excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ concentration of 115 ppm, bringing the combined mean atmospheric (including volcanic) ${\mathrm{CO}}_{\mathrm{2}}$ concentration tree exposure to potentially around $\sim \mathrm{520}$ ppm. For B. nitida this equation yields 133 ppm of mean excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ at the highest flux location, for a combined total mean of potentially $\sim \mathrm{538}$ ppm ${\mathrm{CO}}_{\mathrm{2}}$. 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 ${\mathrm{CO}}_{\mathrm{2}}$ that these trees are exposed to. A ${}^{\mathrm{14}}\mathrm{C}$ tree ring study at Mammoth Mountain found an average yearly volcanic excess ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ than the tree in the Mammoth Mountain study. Additional measurements of tree core ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ and associated soil ${\mathrm{CO}}_{\mathrm{2}}$ fluxes would help corroborate our observations, which were based on a limited number of data points. Though tree ring ${}^{\mathrm{14}}\mathrm{C}$ content in volcanically active areas has been linked to variations in volcanic ${\mathrm{CO}}_{\mathrm{2}}$ emissions, and comparing patterns of ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ to ${}^{\mathrm{14}}\mathrm{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 ${}^{\mathrm{14}}\mathrm{C}$ dating of trees exposed to naturally isotopically distinct excess ${\mathrm{CO}}_{\mathrm{2}}$ is, in fact, unfortunately not a reliable method for these environments due to the well-known ${\mathit{\delta }}^{\mathrm{14}}\mathrm{C}$ deficiency in trees exposed to excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ which is isotopically “dead” with respect to ${}^{\mathrm{14}}\mathrm{C}$, creating spurious patterns that preclude dating by ${}^{\mathrm{14}}\mathrm{C}$ (e.g., Lefevre et al., 2017; Lewicki et al., 2014).

Our data demonstrate that ${\mathrm{CO}}_{\mathrm{2}}$ fluxes through the soil may be a representative relative measure for ${\mathrm{eCO}}_{\mathrm{2}}$ exposure of overlying tree canopies. Forest canopy exposure to volcanic ${\mathrm{CO}}_{\mathrm{2}}$ will vary over time, as will volcanic ${\mathrm{eCO}}_{\mathrm{2}}$; once emitted through the soil into the sub-canopy atmosphere, the gas experiences highly variable thermal and wind disturbances which significantly affect dispersion of ${\mathrm{CO}}_{\mathrm{2}}$ on minute-to-minute, diurnal, and seasonal timescales (Staebler and Fitzjarrald, 2004; Thomas, 2011). These processes cause in-canopy measurements of ${\mathrm{CO}}_{\mathrm{2}}$ concentration to be highly variable, making instantaneous concentration measurements in a single field campaign not representative of long-term relative magnitudes of ${\mathrm{CO}}_{\mathrm{2}}$ exposure. Soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ into the sub-canopy atmosphere as soil ${\mathrm{CO}}_{\mathrm{2}}$ fluxes is therefore expected to better represent long-term input and exposure of tree canopies to ${\mathrm{eCO}}_{\mathrm{2}}$ than direct instantaneous measurements of sub-canopy ${\mathrm{CO}}_{\mathrm{2}}$ concentration. Previous studies at Turrialba have shown that local volcanic ${\mathrm{CO}}_{\mathrm{2}}$ flux is relatively constant on monthly to yearly timescales (de Moor et al., 2016). Therefore, current soil ${\mathrm{CO}}_{\mathrm{2}}$ fluxes should give relatively accurate estimates of ${\mathrm{CO}}_{\mathrm{2}}$ exposure over time. This paper corroborates that expectation by demonstrating strong spatial correlations between volcanically enhanced soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ and volcanic ${\mathrm{CO}}_{\mathrm{2}}$ fluxes but focused on the difference between trees killed by extreme ${\mathrm{CO}}_{\mathrm{2}}$ conditions and those that were still alive (Biondi and Fessenden, 1999). They concluded that the changes in ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ that they observed were due to extreme concentrations of ${\mathrm{CO}}_{\mathrm{2}}$ (soil ${\mathrm{CO}}_{\mathrm{2}}$ concentrations of up to 100 %) impairing the functioning of root systems, leading to closure of stomata and water stress (Biondi and Fessenden, 1999). ${\mathrm{CO}}_{\mathrm{2}}$ does not inherently harm trees, but the extreme ${\mathrm{CO}}_{\mathrm{2}}$ concentrations (up to 100 % soil ${\mathrm{CO}}_{\mathrm{2}}$) 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurements and that variations in our ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurements are more likely to be caused by the direct photosynthetic incorporation of isotopically heavy volcanic ${\mathrm{CO}}_{\mathrm{2}}$. Our ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurements have no statistically significant correlation with stomatal conductance, which suggests that our heavier ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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_{\mathrm{v}}/F_{\mathrm{m}}$ $\sim \mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ fluxes by several orders of magnitude (well over 10 000 g m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$) than the areas we sampled (up to 38 g m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$), making it much more likely that stress from soil acidification causes stomatal closure and affects wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values are best explained by direct photosynthetic incorporation of isotopically heavy volcanic ${\mathrm{CO}}_{\mathrm{2}}$. To the best of our knowledge, this is the first time that a direct correlation between volcanic soil ${\mathrm{CO}}_{\mathrm{2}}$ flux and wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{eCO}}_{\mathrm{2}}$

Short-term plant functional responses at the leaf level to elevated ${\mathrm{CO}}_{\mathrm{2}}$ were highly species-dependent. B. nitida had no statistically significant functional responses to soil ${\mathrm{CO}}_{\mathrm{2}}$ flux, and O. xalapensis only had a weak negative correlation between soil ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$, 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$. 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$) 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${}^{-\mathrm{1}}$, 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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ over the lifetimes of the trees sampled. Due to the continuous nature of the volcanic ${\mathrm{CO}}_{\mathrm{2}}$ enhancement, we are not investigating and analyzing transient events, and our results instead represent spatial variability in excess ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ (Suess effect) to be relevant, if it does affect our samples, it would be beneficial for detection of volcanic ${\mathrm{CO}}_{\mathrm{2}}$, as the Suess effect is gradually increasing the gap between atmospheric and volcanic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$. Since our ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values likely represent several years of growth, small-scale temporal variations in excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ measurements, we would expect this to be represented by some correlation between DBH and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$, though we did find a strong correlation between the completely independent diffuse excess (volcanic) ${\mathrm{CO}}_{\mathrm{2}}$ flux and wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ histories and to study its potential fertilization effect, due to the completely independent nature of the volcanic excess ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ emissions approach to study tropical ecosystem responses to ${\mathrm{eCO}}_{\mathrm{2}}$ – 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 ${\mathrm{eCO}}_{\mathrm{2}}$. 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 ${\mathrm{CO}}_{\mathrm{2}}$. Since the excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$-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 ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$-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 ${\mathrm{CO}}_{\mathrm{2}}$ fluxes and ecological disturbances, which may be utilized to analyze the long-term effects of enhanced levels of ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ and volcanic ${\mathrm{CO}}_{\mathrm{2}}$ could be a game changer by establishing long-term histories of volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${}^{\mathrm{14}}\mathrm{C}$ content have been shown to correlate well with variations in volcanic ${\mathrm{CO}}_{\mathrm{2}}$ flux (Evans et al., 2010; Lefevre et al., 2017; Lewicki and Hilley, 2014), ${}^{\mathrm{14}}\mathrm{C}$ is relatively expensive to measure, limiting the spatial and temporal coverage of data that can be acquired. ${}^{\mathrm{13}}\mathrm{C}$ is an inexpensive alternative to ${}^{\mathrm{14}}\mathrm{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 ${}^{\mathrm{13}}\mathrm{C}$ and ${}^{\mathrm{14}}\mathrm{C}$ in plants that have incorporated volcanic ${\mathrm{CO}}_{\mathrm{2}}$, strengthening the potential for using ${}^{\mathrm{13}}\mathrm{C}$ in this type of study. Further development of the ${}^{\mathrm{13}}\mathrm{C}$ approach to tracking volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${}^{\mathrm{14}}\mathrm{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 ${\mathrm{eCO}}_{\mathrm{2}}$, 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 ${\mathrm{CO}}_{\mathrm{2}}$ with long-term changes in ${\mathit{\delta }}^{\mathrm{13}}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$, diffusively cold emitted through soils into overlying forests. These isotopically heavy volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ into biomass. Each tree studied was co-located with a soil ${\mathrm{CO}}_{\mathrm{2}}$ flux measurement, and their soil ${\mathrm{CO}}_{\mathrm{2}}$ flux signals vary spatially around a continuous long-term local natural excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ source, which creates a local ${\mathrm{CO}}_{\mathrm{2}}$ gradient within which all the sampled trees are found. The excess volcanic ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values in wood appear not to have significantly affected our measurements, indicating that the heavier wood ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ values could be caused by photosynthetic incorporation of volcanic excess ${\mathrm{CO}}_{\mathrm{2}}$. One of the three species studied (A. acuminata) has consistent positive correlations between instantaneous plant function measurements and diffuse ${\mathrm{CO}}_{\mathrm{2}}$ flux measurements, indicating that short-term variations in elevated ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ to study ecosystem responses to elevated ${\mathrm{CO}}_{\mathrm{2}}$ and to use trees as sensors of the changing degassing behavior of volcanic flanks, which is indicative of deep magmatic processes.

Data availability

Data can be found in Tables S1 and S2 in the Supplement or can be requested from Florian Schwandner ([email protected]).

The supplement related to this article is available online at: https://doi.org/10.5194/bg-16-1343-2019-supplement.

Author contributions

FMS and JBF designed the study, and RRB, FMS, JBF, and ED conducted the field work and collected all samples and data with some of the equipment borrowed from GN, who helped interpret the results. TSM processed the samples for analysis. JPL conducted the ${\mathrm{SO}}_{\mathrm{2}}$ analysis, wrote the related methods subsection, and helped interpret the results. VY modeled the anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions, wrote the related methods subsection, and helped interpret the results. CAF created the combined figure showing the ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{SO}}_{\mathrm{2}}$ results and assisted in writing the paper. RRB wrote the publication, with contributions from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We are grateful for LI-COR, Inc. (Lincoln, NE, USA) providing us a loaner ${\mathrm{CO}}_{\mathrm{2}}$ sensor for field work in Costa Rica. We thank Rizalina Schwandner for engineering assistance during sensor integration, OVSICORI (Observatorio Vulcanológico y Sismológico de Costa Rica – the Costa Rican volcano monitoring authority) for logistical and permit support, SINAC (Sistema Nacional de Áreas de Conservación – the Costa Rican national parks service) for access to the Turrialba volcano, and Marco Antonio Otárola Rojas (Universidad Nacional de Costa Rica – ICOMVIS) for invaluable help in the field. We also thank three anonymous reviewers, Akira Kagawa, and the handling editor for very helpful and insightful suggestions that led us to improve the paper. Incidental funding is acknowledged from the Sherman W. Hartman Memorial Fund at Occidental College for funding Robert R. Bogue's field expenses as well as the Jet Propulsion Laboratory's YIP (Year-Round Internship Program) and the Jet Propulsion Laboratory Education Office for funding and support for Robert R. Bogue. Florian M. Schwandner's UCLA contribution to this work was supported by Jet Propulsion Laboratory subcontract 1570200. Part of the research described in this paper was carried out at the Jet Propulsion Laboratory at the California Institute of Technology under a contract with the National Aeronautics and Space Administration.

Review statement

This paper was edited by David Gillikin and reviewed by Akira Kagawa and three anonymous referees.