Landsat-Based Indices Reveal Consistent Recovery of Forested Stream Catchments from Acid Deposition

: Central European forests suffered from severe, large-scale atmospheric depositions of sulfur and nitrogen due to coal-based energy production during the 20th century. High deposition of acid compounds distorted soil chemistry and had negative effects on forest physiology and growth. Since 1994, continuous data on atmospheric deposition and stream runoff fluxes have provided evidence of ecosystem recovery from acidification. In this study, we combined for the first time mass budget data (sulfur deposition and total dissolved inorganic nitrogen (DIN) export) from the GEOMON monitoring network of headwater catchments with annual trajectories of vegetation indices derived from Landsat remote sensing observations. Time series of selected vegetation indices was constructed from Landsat 5, 7, and 8 using Google Earth Engine. Linear regression between the field data and vegetation indices was analyzed using R software. Biogeochemical responses of the forested catchment to declining acid deposition (driven by SO 2 emission reduction) were consistent across all catchments covering various forest stands from different regions of the Czech Republic. Significant correlations were found with total sulfur depositions, suggesting that the forests are continuously and consistently prospering from reductions in acid deposition. Disturbance index (DI) was the only vegetation index that was well-related to changes in forest cover associated with salvage loggings (due to the forest decline) during the 1980s and 1990s. A significant relationship (R 2 = 0.82) was found between the change in DI and DIN export in stream water. Regrowth of young forests in these highly affected areas tracks the most pronounced changes in total DIN export, suggesting a prominent role of vegetation in nitrogen retention. With the Landsat-derived DI, we could map decennial changes in forest disturbances beyond the small scale of the catchments to the regional level (demonstrated here for two protected landscape areas). This analysis showed the peak in forest disturbances to have occurred around the mid-1990s, followed by forest recovery and regrowth. Despite the improvement in forest ecosystem functioning over the past three decades in mountainous areas, emerging threats connected to changing climate will shape forest development in the near future.


Introduction
Forests cover about 40% of the European Union (EU) territory [1], are regarded as the largest terrestrial ecosystem type, and accommodate a large proportion of European biodiversity. Forests provide a multitude of ecosystem services to humans, including climate regulation, timber production, water supply and regulation, clean air, and many others [2,3]. During the past three centuries, European forests have gone through fundamental changes in species composition, and most of the tree populations have been strongly affected by human activities [4]. In the Central European region, forests were largely converted to Norway spruce (Picea abies L., Karst) monocultures. Norway spruce is one of the economically important tree species, but, at the same time, it is vulnerable and susceptible to various biotic, abiotic, and anthropogenic disturbances when planted outside its natural habitats [5][6][7]. During the 20th century, Central European forests suffered from the severe, large-scale deposition of sulfur (S) and nitrogen (N) [8][9][10] that distorted soil chemistry and had negative effects on tree physiology and growth [11][12][13][14]. Among the highest deposition loads are those observed in the area of the Czech, Polish, and German borders, known as the "Black Triangle" [15] due to its high concentration of coal combustion power plants and topography that favored prolonged inversion events. Nevertheless, the Czech Republic is part of a region with the most pronounced decrease in regional SO2 emissions [10,16].
In Central Europe, the retreat of acid deposition started already in the mid-1980s and accelerated in the 1990s after the implementation of "clean air" legislation. The decline of regional SO2 emissions by 90% was followed by a consistent regional decline in sulfur (S) deposition [17]. Despite the modest reduction of N deposition in bulk precipitation, throughfall N fluxes remained generally unchanged. Nevertheless, acid-sensitive catchments started to recover from anthropogenic acidification, and nitrogen inorganic losses were substantially reduced [18], following the retreat of acid deposition, pointing towards better ecosystem immobilization of this limiting nutrient [19]-either by forest vegetation and/or by soil microbes.
In order to study the large-scale effects of reduced acid deposition on soil and forest recovery, a long-term observation network of 14 forested stream catchments, known as GEOMON (GEOchemical MONitoring), was established in the Czech Republic in 1994 [20]. Data from the monitoring program have served numerous studies to provide evidence for stream recovery from acidification [18,21,22] and reveal declines in stream eutrophication [23], increases in stream dissolved organic carbon [24], and effect of climate change on runoff patterns [25]. None of the studies to date have linked the long-term field observations with remote sensing (RS) data, which was the primary goal for this study.
Although the term "forest health" is frequently used, its definition is complex and might be ambiguous [26]. Despite this, many satellite-based RS methods for monitoring of forest health indicators have been developed (reviewed by [27,28]). According to Coppin et al. [29], those methods offer (1) consistent and repeatable procedures, (2) less expensive ways of obtaining information about land surface characteristics, and (3) possibilities for incorporating information from nonvisible parts of the electromagnetic spectrum. Landsat is one of the most used satellite missions for studying forest health properties [26], disturbance [30], and recovery [31]. Its main advantages are good coverage (global coverage with a revisit time of 16 days), sufficient spatial resolution (30 m pixel size), and, perhaps most importantly, its rich archive containing all images acquired since 1972 that allows for studying spatiotemporal changes in ecosystems for periods longer than 30 years [32,33]. Data from Landsat in this field have been incorporated successfully into studies dealing with forest health and structure on various scales ranging from global [34] to continental [35,36] and regional [31,37].
Many studies evaluating forest disturbances have been carried out using Landsat time series [38][39][40]. Usually, these take into account such large-scale ecological factors as bark beetle outbreaks [41], forest wildfires [36], windthrows [42], and harvesting [36,43]. Much less information can be found about the effects of disturbances caused by changes in biogeochemistry, such as soil acidification related to air pollution [44]. These continuous and subtler disturbances can be difficult to assess via RS because there are many environmental variables altering the RS product.
For studying forest disturbances, a disturbance index (DI) was proposed by Healey et al. [45]. It uses Tasseled Cap transformations of Landsat bands and their combinations and has been tested in both Russia and the USA, where it produced more accurate change classification than did original Landsat reflectance data. The DI was subsequently used in other studies located in regions with forest disturbances. Using linear regression, Eshleman et al. [30] found a relationship between the change in total dissolved nitrogen concentration and change in forest DI. A slightly different approach was carried out by Deel et al. [46], who introduced a cumulative version of DI for assessing the impacts of disturbance on the functioning of a forest. They found that this method produced significant relationships with percent canopy cover at plot level and with canopy percent N at catchment level.
The objective of the present study was to evaluate trends in Landsat RS-based indicators in relation to long-term field observations of total sulfur deposition (throughfall fluxes), as a regionally and temporary consistent measure of acid deposition and dissolved inorganic nitrogen (DIN) stream export, as a measure of nutrient retention capacity of forested catchments. We hypothesized that variation in stream water DIN leaching should track i) forest disturbance severity and ii) the extent of forest recovery from acidification.

Study Sites
This study focused on small stream catchments distributed within the GEOMON network [18,20], established in 1994 to study changes in precipitation and stream chemistry due to reductions in SO2 and NOx + NH3 emissions and their impacts on forest ecosystem functioning. Currently, data from 23 years of continuous monitoring are available to study ecosystem biogeochemical changes [18]. The network consists of 15 catchments, representing a variety of environmental conditions in the Czech Republic. Two catchments (Litávka (LIT) and Spálenec (SPA)) were excluded from this study due to limited temporal data records. Therefore, 13 catchments were considered. Their names (and three-letter abbreviations) are presented, and characteristics described in Table 1, and their locations are shown in Figure 1. The catchments are predominantly covered by Norway spruce (Picea abies L. Karst) forests. The only catchments with substantial coverage of broadleaf forest (European beech, Fagus sylvatica L.) are JEZ and LES. Minor representation of broadleaf trees can be found in the POM catchment. The catchments' sizes range from 22 ha (PLB) to 262 ha (MOD), and they are situated mainly in montane and submontane areas (often being parts of protected landscape areas). Most of the catchments are climatically and lithologically predisposed to be sensitive to acidification, with the exception that PLB is well buffered by its ultrabasic bedrock [47]. The three catchments UDL, UHL, and JEZ were most exposed to direct air pollution in the second half of the 20th century and suffered from large-scale deforestation that reached its maximum in the mid of 1990s [44,48].

Field Data
Input fluxes (atmospheric deposition) and stream water runoff fluxes (output) were regularly monitored on a monthly basis using standardized field and analytical methodology in the GEOMON network. Regular monitoring at all catchments started since 1994; however, on some sites (JEZ, LYS, PLB, SAL), measurements started even earlier. Both open area (bulk) and throughfall precipitation were collected using polyethylene collectors. Stream water runoffs were computed by means of Vnotch weirs or flumes with established rating curves. Although water samples were analyzed for detailed chemical composition (summarized in Oulehle et al. [18]), for the purposes of this study, we focused solely on total sulfur deposition derived from the throughfall precipitation fluxes (further denoted as S-input), collected bellow tree canopies, and DIN stream export (further denoted DINoutput). Ammonia concentration contributed minor to the inorganic nitrogen fluxes in stream water; therefore, DIN export equaled NO3 -export. Sulfate and nitrate concentrations were measured by ion chromatography (Knauer 2000). Monthly observations were aggregated into annual fluxes (kg ha −1 year −1 ) covering the period 1994-2016, by multiplying the volume-weighted concentration of solutes by water fluxes.

Landsat Data and Annual Trajectories
The processing workflow is summarized in Figure 2. The time series of Landsat images for 13 studied stream catchments was accessed and analyzed using Google Earth Engine (GEE, [49]; https://earthengine.google.com). We used existing GEE image collections of surface reflectance products from Landsat-5 thematic mapper (TM), Landsat-7 enhanced thematic mapper plus (ETM+), and Landsat-8 operational Landsat imager (OLI). Landsat 5 had started acquiring images in 1984, preceding our field observations by 10 years. We used this time period to assess the prior behavior of the forest ecosystems before the field observations within the GEOMON network were begun and to analyze decennial changes in regional forest conditions starting from 1984.
Landsat images were first filtered according to the following criteria: i) images only from the peak vegetation season June to September were included, ii) cloudy or partially cloudy images were removed, iii) images with more than 25% of no data were removed from the study (this applied to Landsat 7 data only because of the scan line corrector failure). The numbers of valid images varied substantially between the catchments (Table 1). We had only 69 images for MOD, which is located in the top parts of the Giant Mountains (Krkonoše) and has frequent cloud cover, but as many as 336 images for ANE, which is a small catchment located in the central part of the Czech Republic having less cloud cover, and additionally is covered by two Landsat footprints (path/row: 191/25 and 191/26). There was a considerable gap in available Landsat images over the Czech Republic between 1996 and 1999 caused by data transfer regulations that affected many regions of the world [50]. All available images per catchment were used to compute zonal statistics for each year, consisting of median and standard deviation for each band and selected vegetation indices (Table 2). Median was used instead of mean because it was not susceptible to outliers, as could be the case when dealing with large amounts of RS data. In this way, annual trajectories for the period 1994-2016 were composed of each catchment from spectral reflectance and vegetation indices that were selected to be sensitive to forest health changes. One of those vegetation indices (VIs) evaluated was the DI proposed by Healey et al. [45]. It was computed from a Tasseled Cap transformation, consisting of a linear combination of six Landsat bands, into the three components-brightness (B), greenness (G), and wetness (W). The DI was calculated as a linear combination of those three components after they had been first rescaled (Br, Gr, Wr) using their mean (Bµ, Gµ, Wµ) and standard deviation (Bσ, Gσ, Wσ) values computed from reference populations of pixels (Equations 1-4).
The reference populations in this step consisted of pixels for all catchments in 1984. This was done so that the DI could be compared throughout all the catchments with the forest condition in the first year of images acquisition. Although on a larger scale, DI should be computed only over forested pixels, in our study, we regarded all pixels as forested because all catchments with the exception of MOD are almost completely covered by tree canopy.
Other vegetation indices (Table 2) were chosen for their sensitivity to different canopy characteristics. These included the widely used normalized difference vegetation index (NDVI), which is sensitive to canopy greenness; enhanced vegetation index (EVI), which was designed to be more sensitive to changes in areas having high biomass and to reduce the influence of atmospheric conditions; plant senescence reflectance index (PSRI), which is sensitive to leaf senescence by estimating the carotenoids and chlorophyll ratio; transformed chlorophyll absorption reflectance index (TCARI), which is sensitive to chlorophyll pigments. Table 2. Summary of vegetation indices computed from Landsat image data in this study (B -blue, G -green, R -red, and NIR -near-infrared band).

Data Analysis
First, annual trajectories of VIs were plotted as normalized z-scores and visually compared with trajectories of the field-measured S-input and DIN-output at catchment level. Second, the Pearson correlation coefficient was computed for each catchment between field measurements and each individual vegetation index time series data (p > 0.05). Third, adopting the method of Eshleman et al. [30], we computed the change in DI and field-measured S-input and DIN-output as the difference between the linear trend line values at the end and beginning of the observation period. Using an ordinary least-squares regression, we related changes in forest conditions as measured by the difference in DI and measured sulfur deposition and stream water DIN. The linear trend of those scatterplots and R 2 , indicating the strength of the given relationship, were calculated. All these statistical analyses were done in R software [55].
Furthermore, we used DI to analyze decennial changes in forest condition since the beginning of Landsat observations at the regional level. For this purpose, we selected two protected landscape areas: Eagle Mountains (Orlické Hory) hory, where the UDL catchment is located, and Slavkov Forest (Slavkovský les), where the PLB and LYS catchments are located. The two protected landscape areas were selected because UDL and PLB were identified (see Results) to be the most contrasting catchments in terms of forest condition changes. Cloud-free images during peak vegetation season were chosen in 1984,1994,2003,2014, and 2019 for Eagle Mountains and in 1985Mountains and in , 1995Mountains and in , 2003Mountains and in , 2015, and 2019 for Slavkov Forest. In just one case, we were unable to find a completely cloud-free scene (for Slavkov Forest from around the year 2013), so we chose the least cloudy image from that year. Unlike how we handled DI for small catchments, here we selected only forested pixels by applying a mask derived from forest index. Forest index introduced by Ye et al. [56] was computed from 2004 (Eagle Mountains) and 2005 (Slavkov forest) Landsat images. The mask was created by applying an empirical threshold of the forest index (3.6 for Eagle Mountains, 4.0 for Slavkov forest). By doing this, we were able to directly compare every pixel and, furthermore, we were able to distinguish pixels that had been deforested in the past (direct pixel comparison would not have been possible if we would have used a different forest mask for each individual image because the deforested pixels would in that case not have been included).

Long-Term Trends in Measured Sulfur Deposition and Stream Water Dissolved Inorganic Nitrogen Export
Standard scores for bulk sulfur deposition (S-input) and stream water nitrogen export (DINoutput) across all GEOMON catchments are shown in Figure 3. It showed a strong decline in S-input during the 1990s. The most pronounced decline of sulfur deposition was observed at the JEZ, UDL, and UHL, which were those most exposed to acid deposition. The maximum deposition rate for JEZ was 90 kg S ha −1 year −1 measured in 1995, and, in 20 years, until 2016, it decreased to 13.6 kg S ha −1 year −1 . Similarly, the maximum deposition for UDL was 84 kg S ha −1 year −1 measured in 1998, and the lowest was 9.9 kg S ha −1 year −1 in 2012. At UHL, the maximum was 72 kg S ha −1 year −1 in 1995, and a minimum was 6.1 kg S ha −1 year −1 in 2011. This declining trend was consistently observed also at all other sites. In less polluted regions of the central part of the Czech Republic, the maximum value of 23.7 kg S ha -1 year -1 was observed at ANE in 1996, and it decreased to 6.8 kg S ha -1 year -1 in 2015.
In general, the stream water DIN-output also exhibited a decreasing trend, but the variability between the years and between catchments was larger than for S-input (as indicated by the wider shaded area in Figure 3b than in 3a). Especially at the UHL, UDL, and MOD, the maximum DINoutput values (maximum 13.2 kg N ha −1 year −1 measured in 1994 at MOD, and a maximum of 23.7 kg N ha −1 year −1 measured in 1997 at UDL) were observed in the late 1990s, which well corresponded to forest deterioration in the top parts of the northern mountain ranges [58,59]. Some of the catchments (such as ANE, LKV, and SAL, located in the Bohemian-Moravian Highlands), however, exhibited almost no trend and variation between years in stream water DIN-output.

Long-Term Trends in Remote Sensing Indicators
Correlation coefficients between field measurements and Landsat VIs are presented in Figure 4 (crosses indicate non-significant correlations (p > 0.05)). The correlation coefficients were calculated for the period since 1994 (when the field measurements began). With the RS data, however, we were able to track forest changes since 1984.
Results showed that more significant correlations were found between S-input and VIs than between DIN-output and VIs. The S-input showed a significant correlation with DI (6 catchments), NDVI (10 catchments), and TCARI (10 catchments). Only DI showed a strong correlation (Pearson's r > 0.75) on three catchments (LYS, UDL, and UHL). Fewer correlations were found for EVI and PSRI across the catchments. The most noteworthy catchments were, on the one hand, MOD and ANE, where almost no correlations with the VIs were observed, and, on the other hand, UDL and UHL, where correlations were observed with almost all VIs.
In general, almost no significant correlations were observed between VIs and DIN-output. The significant correlations were mainly identified for DI (LYS, POM, UDL, and UHL). Based on the results of the correlation analysis, the most promising VIs for tracing forest changes in relation to S-input were NDVI and TCARI. Nevertheless, when inspecting annual trajectories of the individual VIs ( Figure 5), NDVI showed a positive increasing trend for the 1984-2016 observation period for all the catchments, regardless of known deforestation at UDL and UHL catchments in the early 1990s. Figure 5 shows trend lines for standard scores between field measurements and Landsat VIs. For the sake of brevity, only seven catchments were shown here, as they well represented the geographic areas in the Czech Republic and overall variety between the catchments. Vegetation indices were shown from 1984, meaning that they preceded field observations of S-input and DINoutput by 10 years. The overall steeper slope could be seen in S-input, which also explained the stronger correlations for S-input seen in Figure 4.
Among all the VIs, we chose, for several reasons, to use DI for further analysis. For DIN-output, none of the indices were particularly successful, but DI showed a more significant correlation than did any of the other indices. What is more important, only DI was able to trace forest decline and partial deforestation at catchments most affected by acid deposition (namely JEZ, UDL, and UHL) that occurred before the uniform field measurements started ( Figure 5). NDVI and TCARI were not sensitive to these changes. PSRI showed some sensitivity, but, overall, DI was the most sensitive index. Furthermore, by using DI, we were able to directly compare our results with those from other studies [30,44], also using this index for assessing forest disturbances. Focusing on the trend change within the entire observation period (1995-2016), we observed a positive linear relationship between the S-input and DI with R 2 = 0.58 (Figure 6a). A similar relationship was found also between DINoutput and DI, but the regression coefficient was higher at R 2 = 0.82 (Figure 6b).  Decennial changes in DI over the two protected landscape areas are presented in Figure 7 (Eagle Mountains with UDL) and Figure 8 (Slavkov Forest with PLB). In Eagle Mountains, we could observe a severe increase in DI that corresponded to larger forest disturbances located mainly at the mountain ridges, where UDL is located, between the years 1984 and 1994. On the contrary, no dramatic changes in forest disturbances were detected in the region of Slavkov Forest, except likely small-scale logging activities. Increased DI in 2013 in the northwest area was attributed to imperfections in Landsat image data, such as cloud cover. We classified the entire area into three categories and analyzed the proportion of disturbed, no change, and growth/recovery areas (summarized in Figure 9). We could observe the increase of disturbed pixels from 26.2% to 34.4% across the area of Eagle Mountains as a whole and from 27.7% to 94.5% in UDL between 1984 and 1994. In 2004, only 2.3% of pixels belonged to the growth/recovery category. By 2013, the proportion of those pixels had risen to 64.3%. Finally, in 2019, the growth/recovery category occupied 80.7% of the UDL catchment area.

Acidification Retreat across GEOMON Stream Catchments
Precipitation acidity reversal was driven almost exclusively by the reduction in SO2 emission, and it was consistent across all monitored stream catchments. The decline in S deposition of almost 90% over the past 25 years is noteworthy and paralleled the consistent declines in stream DIN export and increase in precipitation pH [18]. Despite reductions in N deposition over this period, the magnitude of decreases in DIN stream concentration was larger, indicating increasing immobilization of N in GEOMON catchments [18]. The temporarily increasing N retention could have been caused by i) greater immobilization in forest biomass and/or ii) greater immobilization in the soil environment. If the former of these was the pivotal mechanism, then individual indices derived from RS data should be able to detect sink in forest biomass, as was demonstrated using DI for UDL and UHL catchments. Nevertheless, improvement to soil through root exudates from tree labile carbon supply might have a direct impact on nutrient immobilization in soil by stimulating soil organic matter processing and microbial N retention. The deterioration of tree growth and rapid improvement after a decline in an acid deposition have been recently revealed by extensive dendrochronological research undertaken in areas affected by acid rain [13,60,61].

Landsat Data Suitability
The opportunity to obtain cloud-free Landsat images for such relatively small study areas was generally good, although the number of images available varied greatly among the sites (Table 1). For some catchments located in mountainous areas with frequent cloud cover (such as MOD in the Giant Mountains), fewer suitable images were found. There was also a considerable gap in available Landsat images over the Czech Republic between 1996 and 1999 caused by data transfer regulations that affected many regions of the world [50]. For assessing forest disturbance on a larger spatial scale, however, and as we demonstrated here for two protected landscape areas, it was difficult to find a completely cloud-free scene for the entire area during the main growing season. Although there are other satellite sensors, such as MODIS (with a temporal frequency of 1-2 days), that allow for seasonal dynamics to be evaluated, data availability since 2000, together with the coarse spatial resolution of 250 m, made the MODIS product unsuitable for assessment in this study. The recently launched Sentinel-2 offers great potential, with enhanced spatial, spectral, and temporal resolutions compared to Landsat. Its potential for long-term trend studies is still to be realized because Sentinel-2's consistent data acquisition began only in 2015. The incompatibility between field measurements and RS image data and frequent gaps throughout the year due to increased probability of cloud cover during the autumn and winter periods made it impossible to study both the trend between years and the seasonal variations as suggested, for instance, by Verbesselt et al. [62].
Computation of annual VI trajectories was fully automated in GEE and in a manner such that we were able to implement several measures to ensure a high degree of consistency between the images, cloud filtering, filtering of the Landsat 7 striping effect, and selection of median value that was less susceptible to outliers than mean would have been. Even all these measures, however, could not eliminate variability between the years. This likely could be attributed to inappropriate cloud mask and, thus, the potential inclusion of areas affected by artifacts of atmospheric corrections. Other reasons could be temporal inconsistency (even though we considered only images from the vegetation season (June-September)) or topographic effects [63]. Landsat data collections in the GEE area were processed to level-2 atmospherically corrected surface reflectance (Tier 1). Although Landsat 5-7 is atmospherically processed using the LEDAPS algorithm [64] and Landsat 8 data using the LaSRC algorithm [65], we did not observe any noticeable offset in yearly median values after 2013, for which Landsat 8 data are available.

Trends in Remote Sensing Indicators
Significant correlations were found between total S deposition and VIs. This was especially the case for NDVI, which combined joint temporal development of S deposition with concurrent aging of the forests and forest management interventions (namely large-scale harvests in selected catchments). The decline in S deposition was observed consistently for all catchments, and so was a consistent increase in NDVI, which was attributed to growth in forest biomass. NDVI trends extended by 10 years of RS observations prior to the field monitoring did not detect any known forest decline in the top parts of the Jizera (UHL) and Eagle (UDL) Mountains. We observed neither expected decrease in NDVI nor visible changes in EVI and TCARI, and this could be explained by the rapid growth of herbaceous understory. The general trend of NDVI nevertheless inversely mirrored the trend in S deposition ( Figure 5), suggesting the likely positive role of declining precipitation acidity on forest growth. This view is further supported by a dendrochronological study done in the Giant Mountains [13]. Contrarily, DI trends exhibited a rapid increase (similarly as observable for the PSRI) around the mid-1990s and then a steady decrease corresponding to forest regrowth ( Figure S1 in the online Supplementary Materials). DI was successfully applied by Kupková et al. [44] and Mišurec et al. [57] to study forest recovery in the Ore Mountains (Krušné hory) within the Czech Republic.
When observing a linear relationship between the change of DI and S-input or DIN-output ( Figure 6), we could infer that the majority of catchments are situated in the middle cluster, with the exception of UDL and UHL, for which there was a significant decrease of DI in the studied time period (the DI was high there in 1995 because of the deforestation, which is documented in Figure 7 for UDL). Therefore, the relationships were driven by factors reliant primarily on forest management decisions (e.g., harvest) rather than on continuous improvement in forest status. In contrast, PLB is experiencing a subtle increase in DI due to recent selective logging activities connected to bark beetle infestation.
Our simple linear model, examining the relationship between change in DI and change in stream water DIN, was consistent with the results of Eshleman et al. [30], but it was created from a longer data time series (22 years instead of the 6 years studied by Eshleman). Although it showed an even stronger relationship, in drawing conclusions about the strength of the relationship, we must admit that this was mainly due to the influence of data for a few catchments situated outside the big cluster of catchments, showing no change in DI and stream water DIN export. The small number and area (<200 ha) of catchments did not allow for examining the effect of the spatial distribution of forest disturbances. This effect was assessed by Cowles et al. [66], who analyzed data from a large catchment with many embedded sub-catchments. It would be interesting to analyze the trends while including more catchments experiencing major positive or negative changes in forest conditions. This is likely to happen soon due to a rapidly spreading bark beetle outbreak that is destroying productive spruce forests in the Czech Republic [67]. Therefore, we could expect that catchments at middle altitudes (below 800 m) would sooner or later be largely deforested due to salvage fellings. Catchments located in the Moravian-Bohemian highlands (ANE, LKV, and SAL) were already been influenced during 2019. Unfortunately, field data on their S-input and N-output were not yet available at the time of this study. As stated above (in 4.1), the DIN export revealed more year-toyear variability compared to stream nitrate concentration due to the variable hydrological conditions. Oulehle et al. [18] showed a consistent decline of stream nitrate concentration across GEOMON catchments. Such a general decline that is steeper than the decline in a deposition is consistent with an overall improvement of N immobilization in forested ecosystems, either in vegetation (consistent with NDVI and TCARI observations) or in soil by higher N microbial immobilization [68].
Large-scale analysis of DI illustrated a similar trend of forest regrowth in the Eagle Mountains protected landscape area after the 1990s, as was observed by Kupková et al. [44] for the eastern part of the Ore Mountains. On those maps (Figures 7 and 8), DI was also sensitive to the spatial resolution of the satellite system at the edges between the studied forest and other land cover types, such as bare soil, and, therefore, it is important to carefully select only forest pixels or pixels that had been forested in the past to compare those between the years. In future years, even better results will be obtained by using satellite systems with better spatial resolutions, like Sentinel-2.
An increase in the proportion of growth/recovery class, seen in Figure 9b, could be associated with stricter pollution control policies and desulfurization of coal combustion power plants. Both of these have made reforestation possible in the most affected areas. By contrast, an increase in disturbed pixels, seen in Figure 9d, could be attributed to sparse logging activities that occurred at our study site.

Conclusions
This study combined more than 20 years of field observations of total sulfur deposition and stream water dissolved inorganic nitrogen export from small forested catchments in the Czech Republic with vegetation indices derived from the Landsat time-series data.
Biogeochemical responses of forested catchments to declining acid deposition (driven by a reduction in SO2 emission) were consistent across catchments. The reduction of DIN export was not exclusively driven by a reduction in N deposition, but higher forest ecosystem retention was detected. Consequently, RS indices of forest cover properties showed marked changes over the past three decades. DI was well-related to changes in forest cover associated with forest management activities (i.e., salvage loggings during the 1980s and 1990s). Regrowth of young forests in these highly affected areas tracked the most pronounced changes in DIN export, thus suggesting the prominent role of vegetation in N retention. The observed temporal changes in DI on particular catchments were mirrored on the larger, landscape-scale, thereby pointing to the extent of forested areas historically affected by acid rain. Contrary to DI, the NDVI and TCARI indices linearly and consistently changed across a majority of catchments, independently of the mean age of forests within individual catchments. This might suggest that the forests are continuously and consistently prospering due to reduced acid deposition.
A combination of independent approaches using mass budget data and remote sensing indices complemented by recent dendrochronological studies (e.g., [13,60,69]) have revealed consistent improvement in forest status, including increased productivity and consequent N ecosystem retention, as a response to a reversal in acid deposition.
Despite the improvement in forest ecosystem functioning over the past three decades in mountainous areas, emerging threats connected to changing weather patterns will shape forest development in the near future. Our study demonstrated that RS is a vital approach for detecting adverse impacts of environmental changes on forest functioning at large (landscape) and small (headwater catchment) scales. Satisfactory detection of change magnitude is dependent on the proper selection of RS indices, with DI being most effective for changes in forest cover, whereas NDVI and TCARI are more suitable for detecting continuous changes associated with forest growth and/or productivity.

Acknowledgments:
We would like to thank the anonymous reviewers whose comments have improved this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.