Volcanic Diffuse Volatile Emissions Tracked by Plant Responses Detectable From Space

Volcanic volatile emissions provide information about volcanic unrest but are difficult to detect with satellites. Volcanic degassing affects plants by elevating local CO2 and H2O concentrations, which may increase photosynthesis. Satellites can detect plant health, or a reaction to photosynthesis, through a Normalized Difference Vegetation Index (NDVI). This can act as a potential proxy for detecting changes in volcanic volatile emissions from space. We tested this method by analyzing 185 Landsat 5 and 8 images of the Tern Lake thermal area (TLTA) in northeast Yellowstone caldera from 1984 to 2022. We compared the NDVI values of the thermal area with those of similar nearby forests that were unaffected by hydrothermal activity to determine how hydrothermal activity impacted the vegetation. From 1984 to 2000, plant health in the TLTA steadily increased relative to the background forests, suggesting that vegetation in the TLTA was fertilized by volcanic CO2 and/or magmatic water. Hydrothermal activity began to stress plants in 2002, and by 2006, large swathes of trees were dying in the hydrothermal area. Throughout most of the 1990s, the least healthy plants were located in the area which became the epicenter of hydrothermal activity in 2000. These findings suggest that plant‐focused measurements are sensitive to minor levels of volcanic unrest which may not be detected by other remote sensing methods, such as infrared temperature measurements. This method could be a safe and effective new tool for detecting changes in volatile emissions in volcanic environments which are dangerous or difficult to access.

successfully measured volcanic CO 2 emissions with the OCO-2 satellite, but could only detect a few of the highest degassing volcanoes in the world.Surface-level CO 2 anomalies need to be strongly anomalous to make a detectable impact on a satellite-based atmospheric column measurement of CO 2 due to the high concentrations of atmospheric CO 2 (Burton et al., 2013;Schwandner et al., 2017).This constraint inherently limits this technique to a small number of highly active volcanoes and highlights the need for an alternative way to measure volcanic CO 2 emissions from space.Other satellites such as the GHGSat Claire can detect CO 2 at higher spatial resolutions (Jervis et al., 2021) but require focused targets for observation, limiting their use for consistent monitoring of the ∼1,400 active subaerial volcanoes worldwide.
Plant-focused measurements may offer a viable indicator for the spaceborne detection of volcanic volatile emissions.Determining how plants respond to rising levels of atmospheric CO 2 is also one of the most pressing questions in climate science and ecology (Cox et al., 2013;Friedlingstein et al., 2014).Measuring how plants respond to volcanic volatile emissions could provide crucial data for this effort as well as an effective way to detect degassing at hundreds of active but poorly monitored volcanoes around the world (Cawse-Nicholson et al., 2018;Friedlingstein et al., 2014).An exploratory study at Turrialba volcano in Costa Rica found a correlation between soil CO 2 fluxes and heavier ẟ 13 C values in wood from nearby trees, which suggests that those trees incorporated significant amounts of volcanic CO 2 (Bogue et al., 2018).Tree ring samples from trees downwind of the volcanic plume from Turrialba had shifts in sulfur and carbon isotopes, which is caused by the incorporation of volcanic sulfur and sulfur-induced leaf function impairment, representing an increase in volcanic degassing (D'Arcy et al., 2019).Major Normalized Difference Vegetation Index (NDVI) anomalies preceded flank eruptions at Mt. Etna in Italy and Nyiragongo in the Democratic Republic of the Congo in the early 2000s, demonstrating that plants respond to and can potentially anticipate the location of fissure eruptions (Houlié et al., 2006).Two follow-up studies found that the trees on Mt.Etna are relatively insensitive to changes in temperature or water availability, suggesting that CO 2 fertilization may have contributed to the NDVI anomalies there (Seiler et al., 2017(Seiler et al., , 2021)).Tree ring samples taken from nearby trees not destroyed in the eruption displayed significant 14 C and ẟ 18 O anomalies in the years leading up to the eruption, suggesting an uptake of volcanic CO 2 and H 2 O (Seiler et al., 2017(Seiler et al., , 2021)).
Water, sulfur dioxide, and hydrogen sulfide are the other primary gas species emitted by volcanic systems.Water is the dominant component of volcanic gas emissions, often making up more than 90% of the total gas released (Fischer, 2008;Shinohara, 2008).Similar to CO 2 , it is difficult to measure with remote sensing methods due to the high concentrations of H 2 O in the atmosphere (Girona et al., 2015).Elevated H 2 O degassing likely preceded eruptions at Sabancaya, Mt.Etna, and Nyiragongo, establishing it as a potential target for volcano monitoring operations (Kern et al., 2017;Seiler et al., 2021).Volcanic H 2 O degassing has typically been measured at fumaroles or volcanic vents; and we are unaware of any studies which focus on diffuse H 2 O degassing.Unlike H 2 O, volcanic SO 2 emissions are relatively easy to measure with satellites, and estimates of SO 2 degassing are now available for 106 volcanoes around the world (Carn et al., 2017).At volcanoes with well-developed hydrothermal systems, large quantities of sulfur are scrubbed by dissolution in the hydrothermal fluids while ascending through the crust, leading to areas with nearly pure CO 2 emissions at the surface (Kern et al., 2022;Symonds et al., 2001).Trees can respond to elevated concentrations of CO 2 in a variety of ways, including increasing photosynthesis and specific leaf area (Ainsworth & Long, 2005;Saban et al., 2019).NDVI is a measurement of chlorophyll abundance in a given area, which is a proxy for photosynthesis and tree health and thus has potential for detecting tree responses to changes in local volatile concentrations (Gamon et al., 1995).Volcanic activity can also negatively affect or kill plants, which decreases the NDVI of the area.Sulfur-bearing gas emissions, extremely high CO 2 concentrations, and/or increases in soil temperature can kill large areas of plants (Farrar et al., 1995;Jenkins et al., 2012).All three of these phenomena are typically associated with hydrothermal activity or shallow magma, whereas CO 2 emissions without accompanying sulfur gases or major soil temperature anomalies are usually indicative of the degassing of deeper pockets of magma (Symonds et al., 2001).
Large-scale ground temperature increases have recently been detected via remote sensing at a number of volcanoes prior to eruptions, including those which did not appear to produce any other early warning signs (Girona et al., 2021).Diffuse degassing is often strongly correlated with heat flow at volcanoes and hydrothermal systems, and these temperature anomalies are potentially associated with increased degassing (Chiodini et al., 2005;Finizola et al., 2006;Revil et al., 2011;Zimmer & Erzinger, 2003).Increased CO 2 emissions were detected prior to some of these eruptions, including the 2014 eruption of Pico de Fogo in Cape Verde and the 2009 eruption of Mt.Redoubt in Alaska (Bull et al., 2011;Global Volcanism Program, 2014), strengthening the link between diffuse degassing, diffuse heat flux, and eruptions.
In this study, we investigate the history of the Tern Lake thermal area (TLTA), a recently recognized active hydrothermal area in Yellowstone caldera in Wyoming, USA.Hydrothermal activity was identified from satellite images which captured visible hydrothermal alteration in the soil and a ∼33,000 m 2 area of dead trees accompanied by temperature anomalies detected by thermal infrared imaging (Vaughan et al., 2020).These anomalies were first detected in satellite images from 2000, and from 2006 to 2015, the hydrothermally altered area expanded significantly (Vaughan et al., 2020).The principal goal of this study is to investigate how the trees in this area responded to changes in hydrothermal activity, and to determine if these responses provide additional information about the volcanic activity that thermal infrared and visible light images did not capture.

Study Site Description
Yellowstone is a large, restless caldera system characterized by a highly active hydrothermal system which emits an estimated 45 kt d −1 of CO 2 , 96% of which emanates from acid-sulfate hydrothermal areas (Werner & Brantley, 2003).Tern Lake is located in the northeast sector of Yellowstone just inside the caldera rim in an area which hosts a number of hydrothermal zones (Christiansen, 2001;Vaughan et al., 2020).Activity at the Tern Lake hydrothermal area was only recently recognized from satellite imagery that showed hydrothermally altered soil, dead trees, and infrared temperature anomalies (Vaughan et al., 2020).An aerial photo from 1994 showed that there was a forest with no apparent signs of hydrothermal activity (Vaughan et al., 2020).The hydrothermal area slowly expanded from 2000 until 2005, then started expanding more rapidly for the next 10 years before stabilizing at ∼33,300 m 2 around 2015 (Vaughan et al., 2020).This expansion was primarily identified by measuring the size of the tree kill zone, which is likely caused by some combination of excessive heat in the soil, sulfur gasses, or excessively high CO 2 concentrations in the soil (Farrar et al., 1995;Jenkins et al., 2012;Vaughan et al., 2020).Infrared camera images showed an arc-shaped zone of high temperatures above 50°C; but outside this zone, the soil temperatures generally return to ambient levels within a few meters (Vaughan et al., 2020).Vaughan et al. (2020) characterized this as an acid-sulfate area, a type of hydrothermal zone typically associated with high CO 2 and sulfur gas emissions (Vaughan et al., 2020;Werner & Brantley, 2003).These inferred gas emissions were confirmed by the presence of steaming fumaroles and the smell of H 2 S (Vaughan et al., 2020).The trees present in the hydrothermal zone and the area surrounding it have been exposed to significant quantities of volcanic gases in the last two decades, and Vaughan et al. (2020) recommended radiocarbon dendrochronology to determine the area's history of CO 2 emissions (Cook et al., 2001;Evans et al., 2010;Lewicki & Hilley, 2014;Vaughan et al., 2020).The TLTA is an ideal study site due to the decadal-scale growth of the affected area and the presence of undisturbed and otherwise virtually identical forests in close proximity.The progressive evolution of the area allows us to study plant responses to different types and intensities of volcanic activity, and the homogeneity of the species in the forest reduces the number of factors which could complicate interpretation.

Remote Sensing Methods
The NDVI is one of the most widely used vegetation indices for remote sensing.Chlorophyll reflects near-infrared (NIR) light and absorbs red light; therefore, areas which have high concentrations of chlorophyll will have higher NIR values and lower red values (Richardson, 1977).Healthy plants have higher chlorophyll concentrations than stressed plants, so measuring chlorophyll with NDVI provides information on the general health of plants in an area.NDVI is calculated using the red and NIR bands: This calculation yields unitless values which range from −1 to 1, with higher values indicating healthier vegetation.Dense green forest commonly has values between 0.6 and 1, whereas less healthy vegetation will have lower values.Bare rock typically has values around 0, and water will usually have negative values (Jones & Vaughan, 2010).
We used publicly available Landsat 5 and Landsat 8 image databases provided by the United States Geological Survey (USGS).Each image has a 30 m pixel size and was processed for surface reflectance by the USGS (Vermote et al., 2016).We removed any images which had clouds present in the study area and calculated NDVI values for the remaining 886 images.We also discarded any images where the first background site had an average NDVI value of less than 0.4, as that was the lowest value we observed during any of the primary growing season months (June-September) and it is the closest background area to the thermal area.We assume that any images which have a value lower than this are being affected by snow, wildfires, or other factors which could reduce differences between the background area and study site and make comparisons more prone to outliers and misinterpretation.Our methodology is designed to eliminate as many external factors that could be affecting the vegetation as possible so that we can focus exclusively on the effects of the hydrothermal activity.These constraints removed 701 images from the data set, leaving a total of 185 usable images (see Supporting Information S1 for full list of used images).We used ChatGPT to assist with writing and troubleshooting some of the Python code used to download the images and produce Figures 2-5.We reviewed all code taken from ChatGPT to ensure it was accurate and functioned as intended.
We chose the forest immediately north of the TLTA as the first background area for our study (Figure 1).The second and third background areas are in another similar forest, which is 2 km away from the thermal area.We drew the polygons for the background areas with the intent to avoid any clearings in the forest and to include as much contiguous forest as possible, without including any pixels that were near the edge, as these pixels may contain less dense tree cover or bare ground.The three background areas have areas of 41,700, 80,463, and 73,501 m 2 .The TLTA and the first background site are less than 100 m apart and, prior to the expansion of the hydrothermal area in the early 2000s, were both part of one contiguous and relatively homogeneous forest.Diffuse volcanic carbon dioxide emissions enhance air concentrations of CO 2 within a small spatial area, typically returning to background values within a few tens of meters (Fernández et al., 1998;Lewicki et al., 2017;Marín et al., 2005;Rahilly & Fischer, 2021).Diffuse H 2 O emissions probably dissipate similarly, although we are unaware of any studies which have specifically measured this effect in a volcanic or hydrothermal setting.
Our first background area is more than 60 m from the thermal area, so there is very low probability of any significant elevated CO 2 exposure in the background area.Recent fieldwork by one of the authors (RRB) at the TLTA confirmed that the forests in the area consist almost entirely of lodgepole pine (Pinus cortata), the most common type of tree in Yellowstone National Park (Despain, 1990;Yellowstone Volcano Observatory, 2019).This background site is ideal because it is near the study site, has similar terrain, elevation, slope, and aspect, and both sites had the same tree species prior to the hydrothermal activity.The second and third background sites are ∼2 km east of the TLTA and at similar elevations (∼2,530 m).These sites were chosen due to their proximity to the TLTA and their relatively stable, continuous, and homogenous forest (Figure 1).Neither site has any apparent differences in the typical factors which affect plant health, such as soil type, sunlight exposure, and water availability.The nearest hydrothermal area to either of these background sites is more than 500 m distant, and thus has no effect on the vegetation there.
The differences in NDVI values between the hydrothermal and background areas varied considerably over the time period studied.We define these differences as ΔNDVI: We compared ΔNDVI values of the hydrothermal and background areas using several different statistics for each area, including mean, median, maximum (max), minimum (min), 25th percentile (Q1), and 75th percentile 10.1029/2023GC010938 6 of 12 (Q3).For each of these statistics, we calculated the statistic for each area and then compared the two values using Equation 2. For example, the ΔNDVI mean for the first background area is the mean value of the thermal area subtracted from the mean value of the first background area.When these values are negative, the hydrothermal area has higher NDVI values than the background area, suggesting that the hydrothermal activity is promoting plant health by CO 2 fertilization or increased access to water.When these values are positive, the background area is healthier than the TLTA, suggesting that the vegetation in the TLTA is stressed by hydrothermal sulfur emissions, excessively high CO 2 concentrations, and/or high soil temperatures.Trends in these values are significant too; if the values are positive but decreasing, the thermal area still has lower NDVI values than the background area, but the decreasing difference indicates that plant health in the thermal area is increasing relative to the background area.If the various statistics follow different trends, that suggests that some parts of the hydrothermal area are experiencing different effects than other parts.We assume that in the absence of hydrothermal activity, the hydrothermal area and the background areas should have nearly identical values for most or all of these metrics.

Results
From 1984 to 2000, NDVI values in the background areas and the TLTA typically ranged from 0.40 to 0.70 (Figure 3).Mean NDVI values in the TLTA and background areas decreased slightly in 1988, likely due to the wildfires that swept through Yellowstone that year.The quick rebound in NDVI values during the following year confirms that this area did not burn, but the vegetation was still affected in 1988 by the low precipitation and 10.1029/2023GC010938 7 of 12 potentially also smoke from the fires.Precipitation in Yellowstone was relatively low from 1984 to 1998, averaging 506.1 mm yr −1 (Figure 3) (Menne et al., 2012a(Menne et al., , 2012b)).From 1998 to 2023, average yearly precipitation nearly doubled to 892.8 mm yr −1 .The mean NDVI in the three background areas from 1998 to 2023 was 9.5%-12.7%higher than from 1984 to 1998, likely driven by the increased water availability from precipitation.
All ΔNDVI metrics for the three background areas range from 0.04 to 0.16 at the beginning of the study period, with most metrics near 0.1 (Figure 4).The metrics for each background area from 1985 to 2001 have their own yearly variation, but all three follow the same longer-term trend of decreasing to near 0 or slightly negative values by 2001.After 2001, all metrics start increasing, with ΔNDVI min and ΔNDVI Q1 increasing the most rapidly.
The next major inflection point in the graph occurs in 2006, where the rest of the metrics start increasing much more rapidly and generally continue to increase for the remainder of the study period.The ΔNDVI max was the last metric to begin increasing and increased relatively slowly until 2008.ΔNDVI mean , ΔNDVI median , and ΔNDVI Q3 typically followed nearly identical trends throughout the study period, and until 2003 the ΔNDVI Q1 behaved similarly.The ΔNDVI min and ΔNDVI max mostly followed the same broader trends but with more variation due to the inherent volatility of measuring only a single pixel as opposed to the other metrics which incorporate values from a greater number of pixels.
We also examined the spatial distribution of the pixels containing the highest and lowest NDVI values each year, with a particular focus before 2000, which was the previously recognized start of the expansion of the hydrothermal area.Theoretically, the lowest-value pixels should be trees negatively impacted by hydrothermal activity and could potentially indicate the epicenter of any hydrothermal activity that stressed vegetation.The location of the highest NDVI values could indicate where vegetation was positively impacted by hydrothermal activity, either by CO 2 fertilization or increased access to water.The image from each year which had the highest minimum NDVI value was chosen, which was typically in the middle of July but sometimes as late as the end of September and as early as late June.This approach mitigates any seasonal effects by ensuring that the image chosen is in the peak of the growing season and facilitates comparison between years.For each image, the pixels which fall in the highest and lowest bins in a 10-bin histogram were illustrated (see Figure 5 and Figure S1 in Supporting Information S1).
From 1985 to 1992, the lowest NDVI values were generally in the southernmost part of the hydrothermal area (Figure 5).From 1993 to 2001, the location of the lowest values switched to the northernmost part, except for 1996-1997, where it flipped back to the south.This northern area is the same zone where hydrothermal activity was first detected from satellite thermal infrared images in 2000 (Vaughan et al., 2020).Starting in 1999, one of the lowest value pixels was in the center of the hydrothermal area, which became the main epicenter of hydrothermal activity.This pixel consistently had one of the lowest values from 1999 onwards, consistent with images from Vaughan et al. (2020), which show a bare patch of dead trees emerging in the center of the hydrothermal area starting in 2002.The highest values in the thermal area are typically on the western edge, although they move more to the center and eastern side from 1997 to 1998.In the first background area, the highest values are commonly on the eastern side, while the lowest values move around from year to year with no clear pattern (Figure S2 in Supporting Information S1).

Discussion
We identify three notable periods from 1984 to 2022, which represent different types of activities in the study area : 1984-1999, 1999-2006, and 2006-2022.The TLTA evolved rapidly over the duration of the study period, which provides an opportunity to analyze how NDVI changes in response to different types and intensities of hydrothermal activity.The changes in activity after 2000 are well documented by the visible expansion of the tree kill area and the development of thermal anomalies.We use the NDVI responses to interpret the NDVI variations observed prior to 2000.In the absence of hydrothermal activity, we expect background and thermal areas to have relatively similar values for these metrics.

Moderate Activity: 1984-1999
The first period (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) is characterized by a gradual decline in all ΔNDVI metrics from all background areas.This means that the NDVI of the TLTA was generally increasing relative to the background areas throughout this time period.We interpret this as moderate hydrothermal activity promoting growth throughout the thermal area due to increased access to water or CO 2 fertilization.This interpretation also implies that hydrothermal activity at the TLTA started at least 15 years earlier than previously recognized, and that the thermal anomalies observed in 2000 actually represent a reawakening or intensification of activity at the site, as opposed to the onset of activity.This hypothesis is supported by the occurrence of a series of small earthquakes in the early 1980s within 5 km of the thermal area, which are commonly associated with movement of magmatic fluids (Greenfield et al., 2019;Roman & Cashman, 2006;U.S. Geological Survey, 2022).Some of the later expansion of the hydrothermal area is associated with an earthquake swarm in 2010, so this earlier earthquake swarm may also have signified an increase in hydrothermal activity, which then declined over the following decade (Vaughan et al., 2020).The ΔNDVI metrics for the first background area are relatively stable, but the ΔNDVI metrics for the other two background areas are more volatile.There are major drops in ΔNDVI in 1989ΔNDVI in and 1992ΔNDVI in -1993 for both background 2 and 3, which are much less pronounced in the first background area.The first background area is closer to the thermal area than the other two, so this may indicate that there are some climatic conditions that were more strongly affecting the second and third background areas than the first background and thermal areas.Since all three background areas follow the same trends, with varying degrees of intensity, this could also suggest that there were short-term changes in hydrothermal activity in 1989 and 1992-1993, which promoted plant health in the TLTA.
From 1985 to 1991, the lowest NDVI values were located at the southern tip of the hydrothermal area (Figure 5), indicating that hydrothermal activity during this time was probably centered there instead of the northern and central locations where the post-2000 activity first became visible (Vaughan et al., 2020).The location of the lowest value pixel in the TLTA also started migrating in 1992, moving from the southern tip to the center and then to the northern tip of the area in 1993.Except for 1996-1997, the lowest NDVI values are consistently at the northern edge of the thermal area during this time, which is where one of the first infrared temperature anomalies appears in 2000 (Vaughan et al., 2020).This anomaly was first detected in an image from 2000, whereas the lowest NDVI values began appearing there in 1993, 7 years earlier.This may indicate that vegetation in this area was already beginning to be stressed by hydrothermal activity in the early 1990s, several years before it became apparent in thermal infrared or visible images.Fracture pathways through which hydrothermal gases and fluids travel can seal over time due to mineral deposition (Boudon et al., 1998;Edmonds et al., 2003;Stix et al., 1997).The switch of the epicenter of activity from the south to the north that we observe in 1992-1993 may reflect the sealing of the primary transport pathway that fueled activity in the south.We suggest that vegetation can respond to changes in hydrothermal activity long before it becomes apparent from thermal or visible imagery, highlighting its potential for early detection of changes in volcanic activity.

Intensification of Activity: 1999-2006
The second period is characterized by stronger magmatic CO 2 and/or H 2 O fertilization and the beginning of the expansion of the tree kill area.In 1999, the ΔNDVI metrics for background areas 1 and 3 start decreasing more rapidly and some of the lowest NDVI value pixels start appearing in the center of the TLTA.The location of these pixels becomes the center of the tree kill area in subsequent years and is the part of the TLTA with the highest measured temperatures (Vaughan et al., 2020).This is another example of vegetation responding to hydrothermal activity before it is detectable with other remote sensing methods.By 2001, the ΔNDVI mean metrics for all three background areas dropped below 0, the only year in the study period when this occurred.This is the year after the earliest signs of hydrothermal activity were detected in this area, and it strongly suggests that CO 2 fertilization and/or increased access to magmatic water drove a relative increase in NDVI in the thermal area compared to previous years.After 2001, the ΔNDVI Q1 and ΔNDVI min started increasing for all background areas, reflecting the initiation of major vegetation stress in small parts of the thermal area.The other metrics are less consistent, and most of them stay relatively flat until 2006.This indicates that only a small portion of the thermal area in the center and the northern edge had sufficiently intense heat or sulfur emissions to stress plants during this time, while most of the rest of the thermal area was largely unaffected by these negative impacts.This corresponds well with the visible and satellite imagery, which show that thermal output and the size of the tree kill area increased slowly until 2006.From 2002 to 2008, the lowest NDVI values in the thermal area are located exclusively in the center of the TLTA.This indicates that the epicenter of hydrothermal activity in the TLTA has moved from the northern edge to the middle, where it has remained until the present day.

Activity Peaks: 2006-2022
All ΔNDVI values started rapidly increasing in 2006, reflecting the "rapid growth" period described by Vaughan et al. (2020) and indicating that intensifying hydrothermal activity was severely stressing the vegetation in the thermal area.NDVI values in the thermal area continue decreasing rapidly until 2011, after which they continued decreasing but at a slower rate (Figure 4d).By 2011, the hydrothermal area was nearly barren, so the variations in the ΔNDVI metrics at this point are more likely related to fluctuations in the background areas' NDVI values rather than actual changes in hydrothermal activity.Vaughan et al. (2020) noticed significant new tree growth within the bounds of the tree kill area during their field work in August 2019, indicating that hydrothermal activity may be declining in recent years.The trees present there are still very young and too small to have a major effect on NDVI, but if hydrothermal activity does not re-intensify, we would expect NDVI in the thermal area to continue to rise in coming years and for some of the ΔNDVI metrics to start to decrease as these new trees continue to grow.

Future Directions
This approach holds significant promise for future work, although there are some important issues to consider.This project focused on a relatively small area, but our general approach should be applicable to much larger areas as well.The broad flanks of stratovolcanoes consistently emit volatiles and are commonly forested, making them excellent targets for future projects.Volcanoes with well-developed, relatively homogenous, and healthy forests are the ideal targets.Plants that are under significant stress from normal environmental factors such as droughts, wildfires, or high air temperature are less useful for this approach because negative impacts from volcanic sulfur emissions or high soil temperatures will be difficult to disentangle from other non-volcanic stresses.This approach will be easiest to implement in areas similar to the TLTA, which only have one dominant tree species, as the presence of multiple species can influence NDVI values.The spatial distribution of different tree species could lead to certain pixels having lower or higher NDVI values which are unrelated to volcanic or hydrothermal influences, and this issue must be considered in more complicated forest environments.Most volcanoes in arid climates or at high elevations are unlikely to be viable for this approach due to sparse vegetation, which does limit its utility in some major volcanic regions such as the central Andes.Some volcanoes are also cloudy for large parts of the year, which can make it difficult to acquire sufficient cloud-free images to conduct a robust study.Additionally, the single pixel analysis of the spatial variability of high and low NDVI values highlights the need for higher resolution imagery.
In many cases, the ecological effects of volcanic activity are likely subtle, and therefore easy to overlook for human observers.Major large-scale NDVI anomalies like the ones observed at Mt. Etna and Nyiragongo are easily noticeable without deeper analysis (Houlié et al., 2006), but these two examples probably represent a relatively rare phenomenon.These major NDVI anomalies are most likely caused by increases in soil water sourced from degassing magma (Seiler et al., 2021), but more work needs to be done to explore how the magnitude of these effects may vary under different conditions.Increased CO 2 also increases water use efficiency in plants, which may enhance the positive effects of access to magmatically sourced water (Drake et al., 1997).We suspect that most forested, active volcanoes will have certain areas with subtle effects on the local vegetation, which may require large amounts of data and high-level statistical or machine learning analysis for detection.We suggest that these more subtle effects are the most crucial for future studies, as they hold the most potential for widespread use at volcanoes around the world.There is an extensive body of ecological research focused on how ecosystems react to elevated CO 2 concentrations, which could provide a number of different approaches for expanding our techniques (see reviews by Ainsworth and Long (2005) and Saban et al. (2019)).This study focused on NDVI, a relatively simple metric for understanding vegetation health, but future studies could explore specific effects which have been observed in experimental and natural settings, such as increased biomass, increased tree height, increased specific leaf area, increased leaf starch content, and decreased stomatal conductance (Ainsworth & Long, 2005;Saban et al., 2019).A further goal is to demonstrate the potential for gathering ecosystem-scale data on plant responses to CO 2 to provide new insight into how forest ecosystems respond to rising levels of atmospheric CO 2 .Our understanding of how elevated concentrations of CO 2 affect ecosystems has improved considerably in recent years, but there is still a notable lack of field-based data regarding how ecosystems in non-temperate environments react to these conditions (Walker et al., 2021).Future studies of ecosystem responses to volcanic CO 2 emissions could both benefit from and contribute to this body of knowledge, and we strongly encourage collaboration with ecologists to take advantage of this mutual interest.
Ground-based validation will be crucial for confirming the results of our study.Tree ring studies following the model of Seiler et al. (2021) could confirm our findings and help distinguish whether volcanic CO 2 , H 2 O, or some other previously unconsidered factor is responsible for the trends observed.In situ measurements of vegetation traits and functional responses in degassing areas could help identify the strongest plant responses to different types of volcanic degassing.Future studies could then focus on using satellites which are best suited to detecting these responses.We have refrained from providing quantitative estimates of volatile emissions in this paper, but future projects may be able to develop models with ground truthing data which could make yearly estimates of diffuse volcanic volatile emissions in an area of interest possible.This would require ground-based measurements of diffuse volcanic degassing and nearby vegetation, but could be an effective way to scale up degassing measurements from a small area to an entire volcanic edifice, or to identify areas of degassing to target for future ground-based degassing measurements.The TLTA is a relatively small area, so it should be possible to acquire ground truthing data that would be representative of the entire area.If a future project scaled up this approach to a larger area, such as a large volcanic edifice, then ground truthing select portions of the area of interest could allow for quantification of volatile emissions across the entire edifice.Combining ground-based plant and gas emission measurements with our satellite approach could yield a much more comprehensive understanding of volcanic volatile emissions than either approach individually.

Conclusions
Hydrothermal activity at the TLTA was first recognized from satellite detection of infrared temperature anomalies and visual examination of the expansion of the tree kill zone (Vaughan et al., 2020).Our further investigation of the satellite image record focused on measuring plant responses using NDVI demonstrates that hydrothermal activity in this area may have begun significantly earlier than previously recognized, and that the activity observed in the early 2000s is likely an intensification of ongoing hydrothermal activity.The declining ΔNDVI values in 1984-2001 support the hypothesis that there was hydrothermal activity which was promoting plant health in the TLTA during this period.NDVI is higher in the thermal area than the background areas for all measured statistics in 2001, which is shortly after Vaughan et al. (2020) recognized the start of hydrothermal activity, suggesting that increased hydrothermal activity in 2001 released large quantities of CO 2 and/or H 2 O, which promoted plant health.The later expansion of the hydrothermal activity is well tracked by NDVI, and the intensification of activity in 2006 is reflected by large increases in all ΔNDVI metrics.We demonstrate the value of NDVI as a tool for examining past and current volcanic and hydrothermal activity.We would first like to thank two anonymous reviewers for providing detailed and constructive reviews of this paper.Their suggestions greatly improved the quality of the manuscript, and we appreciate them taking the time to review this paper.We would also like to thank Dr. Marie Edmonds for handling editorial duties for this manuscript.We would like to thank Dr. Florian Schwandner for many productive discussions regarding the interactions between volcanic gases and ecosystems.We would also like to acknowledge the USGS for providing the Landsat 5 and 8 images used in this study.Robert Bogue was supported by a Vanier Canada Graduate Scholarship (#175782) during this project.Dr. John Stix and Dr. Peter Douglas were supported by NSERC Discovery Grants.

Figure 1 .
Figure 1.Location of study site within Yellowstone National Park, Wyoming, USA. Park boundaries (black line) and the caldera rim (orange line) shown in top right.The satellite image shown is a Maxar (Vivid) image captured on 17 October 2022, at a resolution of 0.3 m and was sourced from the ESRI World Imagery basemap.

Figure 2 .
Figure 2. (a) Array of images showing an aerial photo, (b) a satellite photo, and (c, d) Normalized Difference Vegetation Index images of the thermal area and the first background area in 1994 and 2022.(a) Is sourced from the National Aerial Photography Program, (b) is a MAXAR (Vivid) photo from the ESRI World Imagery Basemap, (c) is a Landsat 5 image, and (d) is a Landsat 8 image.

Figure 3 .
Figure 3. Annual precipitation and mean Normalized Difference Vegetation Index (NDVI) values for the three background areas and the thermal area from 1984 to 2022.Dots indicate mean NDVI values for individual dates, and the lines represent yearly averages of these values.The vertical red line indicates the date of the first detected vegetation stress in Vaughan et al. (2020).

Figure 4 .
Figure 4. (a) ΔNDVI metrics over time for the first background area, (b) the second background area, and (c) the third background area relative to the thermal area for the time period 1985-2022.The lines in this figure represent yearly average values.The vertical red line indicates the date of the first detected vegetation stress in Vaughan et al. (2020).Positive values indicate that the vegetation in the background area is healthier than that in the thermal area, and negative values indicate the opposite.See Equation 2 for explanation of ΔNDVI.

Figure 5 .
Figure 5. Colored pixels indicate the location of the pixels with the values in the highest and lowest bins in a 10 bin histogram within the Tern Lake thermal area from selected images from 1991 to 2010.Each pixel is colored based on how many standard deviations the Normalized Difference Vegetation Index (NDVI) value of that pixel is from the mean NDVI value for that image (z-score).These four images were chosen because each is representative of several years of activity.A similar figure which includes 32 images from 1985 to 2018 is in Figure S1 in Supporting Information S1, as well as another figure which shows images from the same dates for the first background area (Figure S2 in Supporting Information S1).