Impact of Volcanic Sulfur Emissions on the Pine Forest of La Palma, Spain

In autumn 2021, the largest volcanic eruption on the island of La Palma in historic records took place. The Canary Islands are of volcanic origin and eruptions have always constituted part of their natural disturbance regime. Until recently, their impacts could not be directly observed and studied. Influence of the emission of phytotoxic gases on biodiversity and ecosystem dynamics was hitherto unknown. The recent eruption is still being intensely monitored. We used Sentinel-2 remote sensing data to analyze the spatial extent and intensity of the impact related to sulfuric emissions, aiming to understand the damage patterns in Canary pine forest. The emissions damaged 10% of that forest and affected 5.3% of the Natura 2000 protected areas. We concluded that this is largely due to the toxic effects of the enormous emissions of SO2. We found a clear correlation between the change in the normalized difference vegetation index (NDVI) and distance from the eruption. This pattern was weakly anisotropic, with stronger damage in southern directions. Counteracting effects, such as ash deposition, were largely excluded by combining NDVI change detection with tree cover density. We expect that vegetation damage will be transient. P. canariensis can resprout after forest fires, where most leaves are lost. Consequently, our assessment can serve as a reference for future ecosystem regeneration.


Introduction
Oceanic islands generally receive great attention in terms of nature conservation. Threats to their ecosystems are mostly linked to human impacts, climate change, and invasive introduced species [1,2]. Natural impacts such as volcanic eruptions, landslides, or earthquakes, are known to be important, often resulting in compound disturbances including cascading effects [3], but are widely ignored due to their low frequency and mainly for the fact that they can rarely be monitored directly.
The legacy of volcanic activity has made the Canary Islands an iconic place for volcanology and geology. The term "caldera" (cauldron), for instance, dates back to the Volcanic eruptions contribute significantly to global SO2 emissions. While they account for 7% of global sulfur, they are the source of 18% of sulfur in the atmospheric column [27]. The eruption of Pinatubo in 1991 reached the stratosphere [28], allowing the long-range spread of aerosols. Volcanic emissions from Africa have been measured above the Amazon rain forest [29], and SO2 from southern Chile reached South Africa [30]. However, due to the mix of anthropogenic, biotic, and volcanic sulfur in the atmosphere, to our knowledge there are few studies on the exact impact of volcanic SO2 on forests. The same is true for persistent volcanic degassing [31].
The spread of SO2 from all sources can be monitored by satellite instruments, such as TROPOMI, which is onboard Sentinel-5P (European Earth Observation program Copernicus) [32]. The TROPOMI measurements indicated total emissions of SO2 assessed to an amount of approximately 1 Tg of SO2 [33]. This ranks the eruption within the top 50 SO2 emissions since satellite measurements yielded quantitative values (1978).
In this study, we scrutinize and trace the impact of sulfuric gases on the unique Canary pine forest. We hypothesize that the huge amount of emitted sulfuric compounds would be reflected in chlorotic changes to the pine canopy and that such impacts may be detected from Earth observations. We furthermore predicted that the predominant wind directions during the individual days of volcanic activity will explain the spatial pattern of vegetation damage related to volcanic gases. As gases are rapidly distributed in different layers of the atmosphere depending on the intensity of volcanic activity, it is evident that the recorded emissions cannot be directly translated into effects on the ground, but rather serve as a rough orientation. However, the amount of emitted sulfur is so high that an influence on terrestrial ecosystems is likely to occur.

Study Area
The island of La Palma is one of the westernmost islands of the Canary archipelago ( Figure 2), with an area of 708 km 2 . Its highest elevation of 2426 m a.s.l. is reached at Roque de los Muchachos. The island is subdivided into an older part with higher elevation in the north, and a younger part in the south. In the northern half of the island, volcanic activity stopped hundreds of thousands of years ago. However, the southern half of the island has The spread of SO 2 from all sources can be monitored by satellite instruments, such as TROPOMI, which is onboard Sentinel-5P (European Earth Observation program Copernicus) [32]. The TROPOMI measurements indicated total emissions of SO 2 assessed to an amount of approximately 1 Tg of SO 2 [33]. This ranks the eruption within the top 50 SO 2 emissions since satellite measurements yielded quantitative values (1978).
In this study, we scrutinize and trace the impact of sulfuric gases on the unique Canary pine forest. We hypothesize that the huge amount of emitted sulfuric compounds would be reflected in chlorotic changes to the pine canopy and that such impacts may be detected from Earth observations. We furthermore predicted that the predominant wind directions during the individual days of volcanic activity will explain the spatial pattern of vegetation damage related to volcanic gases. As gases are rapidly distributed in different layers of the atmosphere depending on the intensity of volcanic activity, it is evident that the recorded emissions cannot be directly translated into effects on the ground, but rather serve as a rough orientation. However, the amount of emitted sulfur is so high that an influence on terrestrial ecosystems is likely to occur.

Study Area
The island of La Palma is one of the westernmost islands of the Canary archipelago ( Figure 2), with an area of 708 km 2 . Its highest elevation of 2426 m a.s.l. is reached at Roque de los Muchachos. The island is subdivided into an older part with higher elevation in the north, and a younger part in the south. In the northern half of the island, volcanic activity stopped hundreds of thousands of years ago. However, the southern half of the island has continued to be active until today. A total of eight eruptions have taken place in historic times and more than 80 volcanic cones have been perfectly identified. Here, the current volcanic activity also took place. The volcanic eruption on the western slope of Cumbre Vieja (at a location called Hoya de Tajogaite or Cabeza de Vaca) on La Palma began on 19 September and lasted until 13 December 2021. With this duration, the eruption became the longest on this island, based on historical chronicles dating back to the 15th century. continued to be active until today. A total of eight eruptions have taken place in historic times and more than 80 volcanic cones have been perfectly identified. Here, the current volcanic activity also took place. The volcanic eruption on the western slope of Cumbre Vieja (at a location called Hoya de Tajogaite or Cabeza de Vaca) on La Palma began on 19 September and lasted until 13 December 2021. With this duration, the eruption became the longest on this island, based on historical chronicles dating back to the 15th century. The climate of La Palma is generally subtropical and oceanic with warm temperatures and less pronounced seasonality. Besides the climatic gradients along elevation, the most pronounced differences are evidently related to elevation and aspect [34,35]. The northeast-facing slopes are exposed to trade winds and are regularly covered by clouds condensing on these high mountains. In contrast, the southwest-facing slopes, such as the one with this volcano, are characterized by low cloudiness, arid to semi-arid conditions, and high temperatures, particularly on the bare ground of young volcanic rocks and lapilli (tephra) fields.
La Palma is known as "La Isla Bonita" or "La Isla Verde", which addresses the lush vegetation that covers its steep slopes due to the favorable moist climate. As accessibility is limited, and human population density is low, a large extent of the island is still covered by natural and semi-natural vegetation, e.g., laurel forest and Morello-Erica thickets, particularly in the isolated valleys along the northern coast, and legume shrub vegetation in the alpine zone. Native trees cannot thrive only in the alpine zone above the treeline and in the most arid parts at lower elevations. As human influences related to land use are minor, most of the woodland is mainly characterized by natural disturbance regimes [36,37].
Due to the geographic isolation of individual islands of the Canary archipelago and of their respective ecosystems, all islands host rare, endemic, and endangered plant species [38]. Canary pine forest on La Palma covers 250 km 2 , of which 208 km 2 (83%) is within The climate of La Palma is generally subtropical and oceanic with warm temperatures and less pronounced seasonality. Besides the climatic gradients along elevation, the most pronounced differences are evidently related to elevation and aspect [34,35]. The northeast-facing slopes are exposed to trade winds and are regularly covered by clouds condensing on these high mountains. In contrast, the southwest-facing slopes, such as the one with this volcano, are characterized by low cloudiness, arid to semi-arid conditions, and high temperatures, particularly on the bare ground of young volcanic rocks and lapilli (tephra) fields.
La Palma is known as "La Isla Bonita" or "La Isla Verde", which addresses the lush vegetation that covers its steep slopes due to the favorable moist climate. As accessibility is limited, and human population density is low, a large extent of the island is still covered by natural and semi-natural vegetation, e.g., laurel forest and Morello-Erica thickets, particularly in the isolated valleys along the northern coast, and legume shrub vegetation in the alpine zone. Native trees cannot thrive only in the alpine zone above the treeline and in the most arid parts at lower elevations. As human influences related to land use are minor, most of the woodland is mainly characterized by natural disturbance regimes [36,37].
Due to the geographic isolation of individual islands of the Canary archipelago and of their respective ecosystems, all islands host rare, endemic, and endangered plant species [38]. Canary pine forest on La Palma covers 250 km 2 , of which 208 km 2 (83%) is within protected areas with restrictions for management and land use. A total of 238 plant species have been recorded for this ecosystem, among which 127 plant species are endemic to the Canaries and the other 44 to only La Palma. Most plant species in this ecosystem (90%) are native. Only a few invasive species are established in this ecosystem characterized by extreme natural disturbances [39]. This high degree of endemicity and naturalness in terms of human impact and diverse species assemblages attributes a high value to the Canary pine forest [40] The Canary pine forest is a major ecosystem in the archipelago [40], which has provided a wide range of ecosystem services, including slope stabilization, timber and pitch production, shedding of needle-litter usable as packing material, maintenance of ground freshwater by cloud condensation (horizontal precipitation), and an attraction for tourists and leisure activities. On La Palma, managed plantations are very scarce. Most of its pine forests are dominated by the endemic species Pinus canariensis C. Sm., showing natural conifer dynamics. Large parts of these ecosystems are protected as Natura 2000 nature reserves. Canary pine is the largest Old World pine species, reaching up to 60 m in height and more than 2.5 m in diameter at breast height (DBH) [41]. However, such impressive trees are rare due to the massive historical felling and exploitation of pine forests [42] Canary pine forests are frequently exposed to severe disturbances such as large wildfires, both natural and anthropogenic in origin [43]. Besides the direct impact of fires, nutrient availability in Canary pine forest soils is also affected [44]. However, P. canariensis is highly adapted to fire impact [45][46][47][48]; it can recover after fire by resprouting from epicormic and root collar buds, shaping a peculiar morphology of fire-impacted trees. The trees survive even intense crown fires. Within a few years, the previous canopy cover can be reached again [49]. The species additionally has very thick bark and a high concentration of carbohydrates in the sapwood parenchyma cells [47,50]. The extremely long needles of this pine aid substantially in freshwater provision [51], by contributing to condensation during the frequent fog (clouds touching the ground). In addition, this ecosystem is naturally rather species-poor. The forest structure and canopy are dominated by P. canariensis, while the understory is composed of herb and shrub species. In consequence, damage to the pine tree canopy has a negative influence on the entire ecosystem. In addition, such damage can be seen as an indicator of pollution impact, which may also affect other species that cannot be detected by remote sensing information.

Data Collection
To assess chlorotic damage to the Canary pine forest after the recent eruption in 2021, we chose Sentinel-2 satellite images, which are well-suited for detecting different types of forest damage, e.g., fire damage [37,52], drought [53], bark beetle infestations [54], and defoliation [55]. To compare forest conditions before and after the eruption, all Sentinel-2 images between 19 September 2021 and 14 December 2021 were processed. Overall, 30 images from Sentinel-2 military grid reference system (MGRS) tiles T27RYM and T28RBS were evaluated. However, both smoke from the eruption and cloud formation due to aerosols left very few images clear enough to be usable. Therefore, two images were selected for a before-and-after comparison: 30 November 2019, a near-cloud-free image representing conditions before the eruption, and an image from 14 December 2021, which was the first Sentinel-2 image available with cloud cover low enough to analyze chlorotic damage to the forest. All images were downloaded at processing Level 2A and pre-processed with the R package "sen2r" version 1.5 [56]. All pixels with the levels "No data", "Saturated and defective", and "Cloud (high and medium probability)" in the included scene classification layer (SCL) were masked out, with a 200 m buffer around the clouds masked out as well.
We also attempted to mask other pixel values, such as those containing thin cirrus clouds. The high aerosol concentration in the atmosphere resulting from volcanic activity ( Figure 3) led to us masking most parts of the island. While these aerosols negatively affect image quality if not masked out, the observed chlorosis of the pine forest was severe enough to render these introduced errors negligible. Masking only the mentioned high to medium cloud probability with a 200 m buffer was therefore chosen as a trade-off between data availability and quality. image quality if not masked out, the observed chlorosis of the pine forest was severe enough to render these introduced errors negligible. Masking only the mentioned high to medium cloud probability with a 200 m buffer was therefore chosen as a trade-off between data availability and quality. The normalized difference vegetation index (NDVI) was calculated for all images. Distance from the eruption craters was calculated from the Copernicus EMS data [13] from 9 December. The position of each forest pixel in relation to the eruption was also calculated from this dataset.
The NDVI was calculated for both scenes (30 November 2019 and 14 December 2021) using the following Equation (1) An NDVI difference map (from here denoted "ΔNDVI") was created by subtracting the NDVI before the volcanic eruption from the NDVI afterwards. We used the ΔNDVI as a proxy for vegetation damage. The ΔNDVI was extracted using the packages "rgdal" 1.5-23 [57] and "raster" version 3.4-5 [58] in R version 3.6.3 [59].
An expert-based land cover classification provided by the Cabildo Insular de La Palma (island council) was used as a forest mask for the Canary pine forest. Since the official land cover classification contains some forest gaps and areas with very low forest cover, which negatively influence the pixel-based analysis, it was combined with the Copernicus tree cover density product [60]. A threshold of 50% minimum tree cover was selected, based on visual inspection of problematic areas. We used a combination of both products to extract the pine forest from the Sentinel-2 images. These natural forests exhibit mono-dominant stands of P. canariensis in the canopy. For this reason, no interference with other tree species needs to be considered. Ash deposits also influence vegetation changes detected with the NDVI. Bare ground of volcanic origin in the study area covered with fresh ash showed lower NDVI values than before. This also influenced the decline in NDVI detected in forest pixels with little canopy cover. The combined forest mask starkly reduced this problem, as it only kept forest pixels with reasonably high canopy cover. Additionally, visual inspection of the true-color images created from the Sentinel-2 scenes The normalized difference vegetation index (NDVI) was calculated for all images. Distance from the eruption craters was calculated from the Copernicus EMS data [13] from 9 December. The position of each forest pixel in relation to the eruption was also calculated from this dataset.
The NDVI was calculated for both scenes (30 November 2019 and 14 December 2021) using the following Equation (1): An NDVI difference map (from here denoted "∆NDVI") was created by subtracting the NDVI before the volcanic eruption from the NDVI afterwards. We used the ∆NDVI as a proxy for vegetation damage. The ∆NDVI was extracted using the packages "rgdal" 1.5-23 [57] and "raster" version 3.4-5 [58] in R version 3.6.3 [59].
An expert-based land cover classification provided by the Cabildo Insular de La Palma (island council) was used as a forest mask for the Canary pine forest. Since the official land cover classification contains some forest gaps and areas with very low forest cover, which negatively influence the pixel-based analysis, it was combined with the Copernicus tree cover density product [60]. A threshold of 50% minimum tree cover was selected, based on visual inspection of problematic areas. We used a combination of both products to extract the pine forest from the Sentinel-2 images. These natural forests exhibit mono-dominant stands of P. canariensis in the canopy. For this reason, no interference with other tree species needs to be considered. Ash deposits also influence vegetation changes detected with the NDVI. Bare ground of volcanic origin in the study area covered with fresh ash showed lower NDVI values than before. This also influenced the decline in NDVI detected in forest pixels with little canopy cover. The combined forest mask starkly reduced this problem, as it only kept forest pixels with reasonably high canopy cover. Additionally, visual inspection of the true-color images created from the Sentinel-2 scenes and compared to the ∆NDVI showed a strong decline in NDVI in trees with visible chlorosis and no or little decline in NDVI in nearby healthy trees. Therefore, ash deposits on needles had little effect on the calculated ∆NDVI values in this study.
Data on sulfur emissions were acquired from the Sentinel-5P TROPOMI product focused on atmospheric concentrations ( Figure 4). Sentinel-5P is contributing to the Copernicus Atmospheric Monitoring Services (CAMS). These services inform us about air quality and pollution, with SO 2 as a main data product (see Table 1). The total amount of SO 2 emitted from the eruption per day was retrieved from the Instituto Geológico y Minero de España (IGME) [33] Forests 2022, 13, x FOR PEER REVIEW 7 of 19 and compared to the ΔNDVI showed a strong decline in NDVI in trees with visible chlorosis and no or little decline in NDVI in nearby healthy trees. Therefore, ash deposits on needles had little effect on the calculated ΔNDVI values in this study. Data on sulfur emissions were acquired from the Sentinel-5P TROPOMI product focused on atmospheric concentrations ( Figure 4). Sentinel-5P is contributing to the Copernicus Atmospheric Monitoring Services (CAMS). These services inform us about air quality and pollution, with SO2 as a main data product (see Table 1). The total amount of SO2 emitted from the eruption per day was retrieved from the Instituto Geológico y Minero de España (IGME) [33]  To estimate sulfur deposition in the study area, the weather data during the eruption were acquired from the Cabildo Insular de La Palma weather station at Cumbre Nueva (28°38′49″ N 17°49′30″ W), situated 5.5 km northeast of the eruption. Due to a technical failure of the radio transmitter at the weather station, only data before 22 November 2021 were available. This station measured wind speed, direction, and other meteorological variables every 10 min. Very few data gaps were present, but one day (30 September 2021) was removed from the dataset due to more than two hours of data being missing. Daily SO2 emissions were divided by the number of daily wind observations ( Figure 5). This combination was plotted using the R package "openair" version 2.8-6 [61].  To estimate sulfur deposition in the study area, the weather data during the eruption were acquired from the Cabildo Insular de La Palma weather station at Cumbre Nueva (28 • 38 49 N 17 • 49 30 W), situated 5.5 km northeast of the eruption. Due to a technical failure of the radio transmitter at the weather station, only data before 22 November 2021 were available. This station measured wind speed, direction, and other meteorological variables every 10 min. Very few data gaps were present, but one day (30 September 2021) was removed from the dataset due to more than two hours of data being missing. Daily SO 2 emissions were divided by the number of daily wind observations ( Figure 5). This combination was plotted using the R package "openair" version 2.8-6 [61].

Data Analysis
Connecting the changes in the NDVI in direct comparison before and after the allowed us to quantify the spatial patterns of vegetation damage. Here, we link thes terns of NDVI change with the probability of deposition of sulfur compounds, de respectively from the recorded daily emissions and from the recorded main wind d tions during these days.
Sentinel-2 sensors deliver appropriate remote sensing information due to their tial resolution (up to 10 m) and revisit time of 5 days [62]. However, it was difficult to cloud-free areas during the volcanic eruption because the emitted ashes, smoke, and resulted in dense cloudiness during this period.
ΔNDVI was extracted for all Canary pine forest pixels from the combined f mask. The distance and direction relative to the eruption were also extracted for all f pixels. We then analyzed the relationship of ΔNDVI to the distance from the eruptio

Data Analysis
Connecting the changes in the NDVI in direct comparison before and after the event allowed us to quantify the spatial patterns of vegetation damage. Here, we link these patterns of NDVI change with the probability of deposition of sulfur compounds, derived respectively from the recorded daily emissions and from the recorded main wind directions during these days.
Sentinel-2 sensors deliver appropriate remote sensing information due to their spatial resolution (up to 10 m) and revisit time of 5 days [62]. However, it was difficult to find cloud-free areas during the volcanic eruption because the emitted ashes, smoke, and gases resulted in dense cloudiness during this period.
∆NDVI was extracted for all Canary pine forest pixels from the combined forest mask. The distance and direction relative to the eruption were also extracted for all forest pixels. We then analyzed the relationship of ∆NDVI to the distance from the eruption, as well as the importance of the relative position regarding the eruption. The dataset was split into three parts: north (315 • to 45 • ), east (45 • to 135 • ) and south of the eruption (135 • to 225 • ). In the western direction (225 • to 315 • ), due to the lava stream only two very small tree stands remained and were therefore excluded. Additionally, we split the ∆NDVI for the southern half of the island (south of the mid-point of Cumbre Nueva ridge) into two parts: east and west of the complete ridge of Cumbre Vieja and Nueva. We then tested if there was a significant difference between forest patches on the wind-or leeward sides of Cumbre Vieja/Nueva. We used the R package "raster" [58] to combine the acquired remote sensing data for analysis.
As a measure of the area of damaged forest, the overall area of forest pixels in our forest mask was calculated. We also calculated the area of pixels with a 10% or stronger decline in NDVI and also with a 20% or stronger decline, then cropped it to the area of Cumbre Vieja and Nueva ridge to exclude pixels with declining NDVI from different unrelated causes on the northwestern coast. We used this area to calculate the percentage of Canary pine forest damaged by sulfur dioxide emissions.
To analyze the relationship between the distance from the eruption and ∆NDVI, generalized linear models (GLMs) were fit to the data. Several data transformations (linear, quadratic, and log) were tried out. We selected the best transformation based on AIC and general model fit, and therefore selected GLMs with the formula ∆NDVI~distance + log(distance). Function "pseudoR2" from the package "DescTools" version 0.99.40 [63] was used to calculate Nagelkerke's pseudo-R 2 value for all GLMs.

Results
A clear spatial pattern emerged from Earth observation of NDVI change when comparing the greenness of the pine forest before and just after the volcanic activity. A large proportion of natural Canary pine forest stands on the western slope of Cumbre Vieja showed chlorotic responses to the eruption. True-color images created from Sentinel-2 revealed that the discoloration of the canopy, which can be seen in the field (Figure 6), is also clearly visible in Sentinel-2 images and can be detected for many parts of the Canary pine forest in the surroundings of the eruption.
Cumbre Vieja/Nueva. We used the R package "raster" [58] to combine the acquired remote sensing data for analysis.
As a measure of the area of damaged forest, the overall area of forest pixels in our forest mask was calculated. We also calculated the area of pixels with a 10% or stronger decline in NDVI and also with a 20% or stronger decline, then cropped it to the area of Cumbre Vieja and Nueva ridge to exclude pixels with declining NDVI from different unrelated causes on the northwestern coast. We used this area to calculate the percentage of Canary pine forest damaged by sulfur dioxide emissions.
To analyze the relationship between the distance from the eruption and ΔNDVI, generalized linear models (GLMs) were fit to the data. Several data transformations (linear, quadratic, and log) were tried out. We selected the best transformation based on AIC and general model fit, and therefore selected GLMs with the formula ΔNDVI ~ distance + log(distance). Function "pseudoR2" from the package "DescTools" version 0.99.40 [63] was used to calculate Nagelkerke's pseudo-R² value for all GLMs.

Results
A clear spatial pattern emerged from Earth observation of NDVI change when comparing the greenness of the pine forest before and just after the volcanic activity. A large proportion of natural Canary pine forest stands on the western slope of Cumbre Vieja showed chlorotic responses to the eruption. True-color images created from Sentinel-2 revealed that the discoloration of the canopy, which can be seen in the field (Figure 6), is also clearly visible in Sentinel-2 images and can be detected for many parts of the Canary pine forest in the surroundings of the eruption. ∆NDVI revealed a clear spatial pattern of the damage to the forest canopy-most trees on the western slope of Cumbre Vieja exhibited a strong decline in greenness. This decline is visible for at least 7 km south of the eruption crater, while forest patches on the outer slopes of the Caldera de Taburiente and Cumbre Vieja to the north of the eruption craters exhibited a weaker decline in NDVI (Figure 7). The decrease in greenness in the vicinity of the eruption crater showed a clear relation to distance; forest patches closer to the crater had a stronger decline in NDVI compared to those further away (Figure 7). Overall, the NDVI decrease due to visible chlorosis was clearly located mostly in the area around the eruption and fits visual field observations of such damage. Pine forests in other areas of the island exhibit similar NDVI values in 2021 when compared to 2019. The few small spots of visible NDVI decline elsewhere on the island cannot be linked to visible chlorosis and most likely had different sources, e.g., atmospheric effects, declining vegetation health for other reasons or image artifacts. Of the 183 km 2 of forest in our forest mask, 18.8 km 2 (10.3%) showed an NDVI decline of 10% or more, while 13.4 km 2 ( 7.2%) underwent a decline of 20% or more. The damaged area lies almost completely within the area protected under Natura 2000.
outer slopes of the Caldera de Taburiente and Cumbre Vieja to the north of the eruption craters exhibited a weaker decline in NDVI (Figure 7). The decrease in greenness in the vicinity of the eruption crater showed a clear relation to distance; forest patches closer to the crater had a stronger decline in NDVI compared to those further away (Figure 7). Overall, the NDVI decrease due to visible chlorosis was clearly located mostly in the area around the eruption and fits visual field observations of such damage. Pine forests in other areas of the island exhibit similar NDVI values in 2021 when compared to 2019. The few small spots of visible NDVI decline elsewhere on the island cannot be linked to visible chlorosis and most likely had different sources, e.g., atmospheric effects, declining vegetation health for other reasons or image artifacts. Of the 183 km² of forest in our forest mask, 18.8 km² (10.3%) showed an NDVI decline of 10% or more, while 13.4 km² ( 7.2%) underwent a decline of 20% or more. The damaged area lies almost completely within the area protected under Natura 2000. . Some minor parts of the forest are not shown in the NDVI product due to cloud cover. Yellow coloration indicates no change, orange and red a decrease in greenness as an effect of the eruption, and green an increase in greenness. Background map: Hillshade generated from a digital elevation model. Figure 8 shows the relationship between ΔNDVI and the distance from the eruption craters for 100,000 random points within the pine forest. A clear pattern is visible, with the highest decline in NDVI, up to 100%, closest to the crater and smaller declines further away. At roughly 7 km, the majority of forest pixels show only limited or no NDVI change due to the eruption. Nonetheless, there are several groups with a strong (30-50%) decline . Some minor parts of the forest are not shown in the NDVI product due to cloud cover. Yellow coloration indicates no change, orange and red a decrease in greenness as an effect of the eruption, and green an increase in greenness. Background map: Hillshade generated from a digital elevation model. Figure 8 shows the relationship between ∆NDVI and the distance from the eruption craters for 100,000 random points within the pine forest. A clear pattern is visible, with the highest decline in NDVI, up to 100%, closest to the crater and smaller declines further away. At roughly 7 km, the majority of forest pixels show only limited or no NDVI change due to the eruption. Nonetheless, there are several groups with a strong (30-50%) decline in NDVI at larger distances, e.g., at 8 and 10 km. A GLM with the formula ∆NDVI~distance + log(distance) showed a significant relationship between ∆NDVI and distance, and had a Nagelkerke's pseudo-R 2 of 0.41.
Forests 2022, 13, x FOR PEER REVIEW 11 of 19 in NDVI at larger distances, e.g., at 8 and 10 km. A GLM with the formula ΔNDVI ~ distance + log(distance) showed a significant relationship between ΔNDVI and distance, and had a Nagelkerke's pseudo-R² of 0.41. The three subsets representing forests north, east, and south of the eruption show only slightly different relationships between NDVI decline and distance from the eruption (Figure 9). Forests south of the eruption had more affected forest patches closer to the crater compared to the areas northward, at least for the nearest 5 km. This is expressed by the fitted curve for the southward direction, which runs consistently below the curves for northward and eastward for distances between 2 and 5 km, indicating more severely damaged forest patches for those distances. The eastern dataset shows a steeper increase, with most forest pixels roughly 4 km away from the eruption already exhibiting only little to no decline in NDVI. The corresponding GLMs are all highly significant and have pseudo-R² values of 0.54, 0.61, and 0.48, respectively. The three subsets representing forests north, east, and south of the eruption show only slightly different relationships between NDVI decline and distance from the eruption (Figure 9). Forests south of the eruption had more affected forest patches closer to the crater compared to the areas northward, at least for the nearest 5 km. This is expressed by the fitted curve for the southward direction, which runs consistently below the curves for northward and eastward for distances between 2 and 5 km, indicating more severely damaged forest patches for those distances. The eastern dataset shows a steeper increase, with most forest pixels roughly 4 km away from the eruption already exhibiting only little to no decline in NDVI. The corresponding GLMs are all highly significant and have pseudo-R 2 values of 0.54, 0.61, and 0.48, respectively. Copernicus EMS from 9 December 2021). The 100,000 random data points extracted from both rasters were split into three subsets: (a) north of the eruption (225° to 315°), (b) east of the eruption (135° to 225°), and (c) south of the eruption (45° to 135°). Apart from two very small tree stands situated west of the eruption, most trees in this direction were covered by lava. This direction was, therefore excluded from the analysis. The red lines represent a hghly significant (***) GLMs of the formula ΔNDVI ~ distance + log(distance) with Nagelkerke's pseudo-R² values of 0.54, 0.61, and 0.48, respectively.

Discussion
The strong impact of the volcanic eruption on the Canary pine forest is difficult to explain by heat, even if fluid lava was hotter than 1000 °C. Heat impact would only be effective in close proximity to the lava. Other impacts related to the volcanic activity, such as emission of CO2 and volcanic ashes, are also unlikely to result in chlorotic discoloration of pine needles. Ash initially covered a large amount of pine trees, but was removed by wind and rain. The evergreen pine forest did not show substantial seasonal variation in spectral traits of the canopy. Seasonal patterns as a response to precipitation, however, can only be expected for the herbaceous vegetation in the understory. As a consequence, canopy density needs consideration to avoid a strong influence from the phenology of field layer vegetation. This was performed in this study, based on a digital tree cover product, we considered only forest surfaces where tree cover was >50%.
In addition to the distance related decay in plant damage, we also found a weak anisotropic distribution with stronger impact southward from the volcano. Indeed, the western lowlands downslope are mainly covered by land use and settlements. The eastern slopes were protected from the eruption by the Cumbre Vieja ridge and, more so, the then prevailing wind direction. The major wind directions during the eruptive phase were ENE and E, which explains this asymmetric pattern of damage. However, overall, the anisotropic relationship was weak, indicating that most of the SO2 emission reaching the pine

Discussion
The strong impact of the volcanic eruption on the Canary pine forest is difficult to explain by heat, even if fluid lava was hotter than 1000 • C. Heat impact would only be effective in close proximity to the lava. Other impacts related to the volcanic activity, such as emission of CO 2 and volcanic ashes, are also unlikely to result in chlorotic discoloration of pine needles. Ash initially covered a large amount of pine trees, but was removed by wind and rain. The evergreen pine forest did not show substantial seasonal variation in spectral traits of the canopy. Seasonal patterns as a response to precipitation, however, can only be expected for the herbaceous vegetation in the understory. As a consequence, canopy density needs consideration to avoid a strong influence from the phenology of field layer vegetation. This was performed in this study, based on a digital tree cover product, we considered only forest surfaces where tree cover was >50%.
In addition to the distance related decay in plant damage, we also found a weak anisotropic distribution with stronger impact southward from the volcano. Indeed, the western lowlands downslope are mainly covered by land use and settlements. The eastern slopes were protected from the eruption by the Cumbre Vieja ridge and, more so, the then prevailing wind direction. The major wind directions during the eruptive phase were ENE and E, which explains this asymmetric pattern of damage. However, overall, the anisotropic relationship was weak, indicating that most of the SO 2 emission reaching the pine forest was omnidirectional. This is supported by the wind direction data. Especially during the early days of the eruption, the wind blew nearly exclusively from the east and northeast, distributing SO 2 towards the west. On many days, wind speeds were also low, causing volcanic gases to rise into higher atmospheric layers before being spread. The 2021 La Palma volcanic eruption was intensely monitored, both in the field and by Earth observation [64]. Many scientists and institutions, such as the Spanish National Geographic Institute (IGN), Canary Islands Volcanology Institute (Involcan), Spanish National Research Council (CSIC), Univerisity of la Laguna (ULL), University of Las Palmas de Gran Canaria (LPGC), Cabildo Insular de La Palma, and the Canary Islands Government were involved, among others. This intense survey is providing invaluable data resources that will improve our understanding of volcanic impacts. In addition, worldwide attention and the European Copernicus Earth Observation program have yielded Big Data and geospatial information.
Here, we focus only on a small section of this treasure, to point out the direct ecosystem responses during and immediately after the eruption. We are convinced that many other studies will follow.
Through Earth observations, negative impacts on the endemic Canary pine forest ecosystems-among which the volcanic cones emerged-can be precisely documented, as in in this study. Evidence suggests that gaseous emissions are the main driver of these responses. Such natural pollution events took place frequently in the Earth's history, but cannot easily be traced back due to their elusive and ephemeral properties [65,66]. Nevertheless, as shown in this study, the effects of such cases of natural pollution can be substantial. In the case of La Palma, 10.3% of the island's pine forests and 5.2% of its protected areas have been affected, and they host many unique endemic species. Currently, it is too early to say whether the small populations of endemic plant species located in the southern half of the island have been severely affected. However, according to field observations, they have not been driven to extinction yet. Only a wild chickpea Cicer canariense A. Santos & G.P. Lewis and a new species currently named Parolinia aridanae, description pending [67], are potentially threatened plants, but to date their conservation is not compromised [68].
The island of La Palma can be seen as a natural laboratory to study ecological, geomorphological, and biogeochemical processes [69,70]. The high degree of naturalness of ecosystems reduces confounding anthropogenic influences. The volcanic eruption at Hoya de Tajogaite in 2021 was the largest such event within historical times on La Palma. Records of previous eruptions date back to 1450, which are still visible as uniquely distinct lava flows. However, the vast lapilli fields on Cumbre Vieja ridge cannot be attributed directly to individual volcanic cones or events. Furthermore, likely chemical impacts in the past are not reflected, protocolled, or preserved. Volcanic gases are transient, but have a strong impact on vegetation near eruptions, which is surprisingly little studied (e.g., [71][72][73]). In 2021, the amount of sulfuric compounds released was precisely protocolled on a diurnal basis. In combination with recorded wind directions, it allows us to assess the spatial influence of sulfuric gases on the ground. No other hampering impacts, such as fire or drought, occurred in this short period. Thus, it is most likely that the NDVI change was caused by the enormous amounts of sulfuric emissions.
In contrast to the widespread literature dealing with fire impact on Canary pine and on the dynamics of its forest ecosystems (e.g., [37,43,44,48,74], research was lacking until now on the impact of volcanic gases on these trees and forests. However, such events occurred frequently during the existence of the archipelago. The impact of pollutants, including sulfuric compounds, on this pine were hitherto only studied in connection with anthropogenic environmental pollution related to car traffic [75]. Exposure to such types of pollution, however, is many orders of magnitude lower than the toxic compounds released from this recent eruption. The total amount of SO 2 emissions during the active phase of the volcano was enormous, even if the real sulfur emissions cannot be measured precisely while released during that phase. Estimates range up to 1 Tg (1,000,000 tons) of SO 2 [33], which makes this event important at a global level.
Atmospheric measurements from space deliver good proxies and estimates for the order of magnitude. Local measurements of air pollution in settlements that aim at health care are less representative of landscape-scale pollution impacts. In this study, we decided to use data from satellite sensors due to their spatial cover, methodological standardization, and daily availability. Data from different sensors deviate due to technological factors (e.g., Sentinel-5P TROPOMI vs. Nasa Suomi-NPP/OMPS). In this study, this is less relevant as we are interested in the spatial patterns of the impact on the ground. A substantial proportion of the emitted sulfur compounds are projected into the high strata of the atmosphere and transported over long distances, without contact with land surfaces and ecosystems. Nevertheless, satellite sensors can inform on the quantitative development of emissions for every single day. This helps to understand the role of prevalent wind directions and to identify processes leading to the anisotropic distribution of pine damage on La Palma, in addition to the mere distance from the source of emission. On a global scale, sensors from the Copernicus program can be combined to monitor the distribution of emitted SO 2 and to investigate forest responses to such events.
The general importance of volcanic eruptions and gas input in atmospheric chemistry has been highlighted [76,77]. Particularly, emissions of CO 2 and SO 2 are seen as substantial contributions to changes in atmospheric chemistry. However, there is little knowledge about the effects of such gaseous emissions on terrestrial ecosystems. The chlorotic and necrotic responses of pine needles are documented as related to atmospheric sulfur exposure [78], but the ecological literature on atmospheric sulfur impacts is mostly focused on industrial pollution and acid rain (e.g., [79]). As a consequence, these studies are strongly biased towards northern-hemisphere conifer forests. However, it is crucial to better understand the role of such temporarily limited events in ecosystem and biodiversity dynamics [66].
Andres and Kasgnoc [80] compiled the first global inventory of subaerial volcanic sulfur emissions, based on SO 2 measurements en situ. Using data from the second half of the 20th century, they estimated an annual average of 13 Tg SO 2 emissions from volcanoes. This worldwide quantity illustrates the importance of the 2021 eruption on La Palma, with approximately 1 Tg SO 2 emissions. Recently, several groups studied volcanic SO 2 emissions with the support of remote sensing [81][82][83] and modeling [84]. However, there are few studies on the impact of volcanic gases on vegetation [31]. The recent volcanic eruption on La Palma provided a chance to study this fundamental process that happened so frequently in the Earth's history. Highly precise monitoring of the volcano during its active phase, together with Earth observation data, quantification of sulfur emissions from space, and detection of vegetation responses has offered an unprecedented opportunity. The pattern identified will also serve as a reference for future changes in the drivers of ecosystem changes and responses in the biodiversity [85,86].

Conclusions
The eruption on the west slope of Cumbre Vieja (La Palma, Canary Islands, Spain) gained global attention. This event developed into the strongest and longest-lasting eruption in the island's historical records. The emission of ashes and gases had a far-reaching impact, besides the large lava field that now covers former settlements, agricultural fields, and semi-natural vegetation, along with two new volcanic deltas (fajanas) on the coast. Here, we focused on the effects of gaseous emissions, in particular SO 2 and its derivate compounds such as sulfuric acid, on the endemic pine forest of the island. This forest is largely protected through Natura 2000 and is of decisive importance for ecosystem services, such as the provision of freshwater and carbon sequestration.
We successfully used remote sensing-based analyses to find a clear spatial pattern that can be explained by distance decay towards the source of eruption. Chlorotic damage to the pine forest canopy was clearly visible in the Sentinel-2 earth observation data. Due to the lack of large forests to leeward of the eruption, we were not able to prove nor disprove the importance of wind direction and only found a weak anisotropic pattern. Stands of P. canariensis are known to exhibit high resilience towards disturbances, providing a unique study ground for testing an emerging theory on pulse dynamics in ecosystems [66]. This tree species re-sprouts after fire impact and total defoliation. For this reason, we expect this forest ecosystem to recover in the following years. This study also serves as comparison for future monitoring of ecosystem recovery. The specific conditions during the active phase released a large amount of acid aerosols that caused cloud formation in addition to the frequent plume or fan of ashes and smoke from the volcano. Nevertheless, this location on the less humid leeward side of the island will provide many cloud-free days in the future to monitor how the affected pine forests recover.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.