Contrasting 20-year trends in NDVI at two Siberian larch forests with and without multiyear waterlogging-induced disturbances

The fate of a boreal forest may depend on the trend in its normalized difference vegetation index (NDVI), such as whether the NDVI has been increasing significantly over the past few decades. In this study, we analyzed the responses of two Siberian larch forests at Spasskaya Pad and Elgeeii in eastern Siberia to various waterlogging-induced disturbances, using satellite-based NDVI and meteorological data for the 2000–2019 period. The forest at Spasskaya Pad experienced waterlogging (i.e. flooding events caused by abnormal precipitation) during 2005–2008 that damaged canopy-forming larch trees and increased the abundance of water-resistant understory vegetation. By contrast, the forest at Elgeeii did not experience any remarkable disturbance, such as tree dieback or changes in the vegetation community. Significant increasing NDVI trends were found in May and June–August at Elgeeii (p < 0.05), whereas no significant trends were found at Spasskaya Pad (p > 0.05). NDVI anomalies in May and June–August at Elgeeii were significantly associated with precipitation or temperature depending on the season (p < 0.05), whereas no significant relationships were found at Spasskaya Pad (p > 0.05). Thus, the 20 year NDVI trend and NDVI–temperature–precipitation relationship differed between the two larch forests, although no significant trends in temperature or precipitation were observed. These findings indicate that nonsignificant NDVI trends for Siberian larch forests may reflect waterlogging-induced dieback of larch trees, with a concomitant increase in water-resistant understory vegetation.


Introduction
Monitoring changes in vegetation and carbon dynamics in response to environmental changes is key to predicting future responses of Earth environments to climate change (Davidson et al 2006, Luo et al 2016, Sato et al 2020. Broadly observed increasing, and less widely observed decreasing, trends in satellite-derived vegetation indices have been documented in pan-Arctic ecosystems in recent decades (Guay et al 2014, Zhu et al 2016, Wang and Friedl 2019, Myers-Smith et al 2020, so-called greening and browning trends, respectively. The normalized difference vegetation index (NDVI), the most popular vegetation index derived from satellite observations based on surface reflectance of red and near-infrared bands (Guay et al 2014, Zhu et al 2016, Wang and Friedl 2019, Myers-Smith et al 2020, is increasing by 0.005-0.010 yr −1 on average in the pan-Arctic region (Guay et al 2014, Wang andFriedl 2019). Increasing NDVI trends are likely caused by rapidly rising temperatures, changes in precipitation, and a gradually increasing concentration of atmospheric CO 2 (Zhu et al 2016).
However, understanding of the patterns and trends of greening and browning is still lacking (Myers-Smoth et al 2020). Guay et al (2014) and Wang and Friedl (2019) did not find any significant NDVI trend in recent decades for 50% or more of the terrestrial ecosystem area in the pan-Arctic region. Elucidating the ecological features of ecosystems with nonsignificant NDVI trends is important for predicting their future trajectories, in particular to determine whether these ecosystems will have the same outcomes as those exhibiting significant trends.
To examine the ecological features of an ecosystem exhibiting a nonsignificant NDVI trend, we investigate a larch forest at Spasskaya Pad near Yakutsk in eastern Siberia. Larch forest comprises the typical boreal forest vegetation in northern Eurasia, in particular in eastern Siberia (Kolchugina andVinson 1993, Sato et al 2020). Larch trees are coniferous but deciduous, and the floor of larch forests is covered with understory vegetation consisting of mosses, lichens, and shrubs. This understory vegetation makes up a significant proportion of boreal forests and thus contributes substantially to the total ecosystem carbon flux (Baldocchi et al 2000, Iida et al 2009, Ikawa et al 2015. The larch forest at Spasskaya Pad experienced waterlogging (i.e. flood events) owing to abnormal precipitation in the summers of 2005-2008 that greatly damaged larch tree canopies and increased the cover of waterresistant understory vegetation, including grasses and shrubs . In contrast to these significant disturbances and changes in the composition of vegetation, eddy covariance flux observations from the forest for the 2004-2014 period revealed no apparent changes in gross primary production at the whole-ecosystem scale . However, the relative contribution of understory vegetation to total gross primary production gradually increases after waterlogging . Thus, the NDVI trend for the Spasskaya Pad forest in recent decades is expected to be nonsignificant, in particular from 2000 to 2019, although the composition of its vegetation changed significantly after waterlogging in 2005-2008. The NDVI for another larch forest at Elgeeii near Ust-Maya, located 300 km southeast of Spasskaya Pad, is expected to exhibit a significant increasing trend, because waterlogging from the large amounts of precipitation in the summers of 2005-2008 had no remarkable effects on the vegetation (i.e. no larch tree diebacks or changes in the vegetation community; Kotani et al 2014). Kotani et al (2014) found several differences in how carbon is taken up at Spasskaya Pad and Elgeeii, such as a difference between the two forests in the seasonal time course of carbon uptake and less efficient use of water at Spasskaya Pad. Hence, the NDVI time course is expected to differ between the two larch forests. In addition, a comparison between the two forests of the NDVI-meteorology relationship will provide insight into the fates of boreal-forest ecosystems with and without a significant NDVI trend over the past few decades.
In this study, we compared NDVI trends over the 2000-2019 period between the forests at Spasskaya Pad and Elgeeii to evaluate a hypothesis that the NDVI trend for the past two decades differs between these forests with and without multiyear waterlogging-induced disturbances. We also assessed relationships between NDVI, temperature, and precipitation to evaluate another hypothesis that the NDVI-meteorology relationship differs between the two forests. If both hypotheses are true, it is possible that some pan-Arctic forest ecosystems with no significant NDVI trends might have experienced significant disturbances and will likely follow a different trajectory under climate change. For the analyses, we used 20 year NDVI time-series data from Moderate Resolution Imaging Spectroradiometer (MODIS) on Terra. The MODIS NDVI is a major vegetation index that has been analyzed in previous studies of pan-Arctic greening and browning (Guay et al 2014, Zhu et al 2016, Wang and Friedl 2019, Myers-Smith et al 2020. We determined trends in the NDVI using the Mann−Kendall test, as in previous studies of pan-Arctic greening (Guay et al 2014, Zhu et al 2016. NDVI-meteorology relationships were evaluated using path analysis, which is an extension of multiple regression analysis that allows for the examination of a possibly covarying relationship between temperature and precipitation. Path analysis is commonly used in studies of terrestrial carbon flux (Saito et al 2009, Nagano et al 2018, in which intrinsic connections among variables may be pronounced.
The Spasskaya Pad forest, but not the Elgeeii forest, experienced substantial changes in the composition of its understory and overstory (trees that make up the forest canopy) vegetation after waterlogging in 2005-2008. Briefly, the density of living larch trees with a diameter at breast height (DBH) ⩾ 5 cm was 824 trees ha −1 in 1997, but decreased to 660 trees ha −1 in 2018. Birch and willow densities increased from 216 and 1244 trees ha −1 in 1998-1724 and 3400 trees ha −1 in 2018, respectively. Nevertheless, the percentage of birch and willow with DBH ⩾ 5 cm was only 3%, even in 2018. In the Elgeeii larch forest, the density of larch trees with DBH ⩾ 5 cm was 740 trees ha −1 in 2008, and the percentage of birch and willow with DBH ⩾ 5 cm was 14%. Thus, both forests were still designated larch-dominated forests after the waterlogging events. The leaf area index (LAI) values at Spasskaya Pad and Elgeeii in the summer of 2012 were 2.4 and 2.1, respectively . Understory vegetation contributed to 63% and 33% of the total LAI at Spasskaya Pad and Elgeeii, respectively, in summer 2012 .

Satellite NDVI data
NDVI data for the growing season (May-September) for 2000-2019 were derived from MOD09GQ (a MODIS product that includes daily surface reflectance from red and near-infrared bands at 250 m resolution), which is available through Google Earth Engine (https://earthengine.google.com, accessed on 24 May 2021). MODIS data from Google Earth Engine were originally prepared by the NASA EOS-DIS (Earth Observing System Data and Information System) Land Processes Distributed Active Archive Center (LP DAAC) of the USGS Earth Resources Observation and Science (EROS) Center (https://lpdaac.usgs.gov, accessed on 24 May 2021). Using the Google Earth Engine API for Python processing, we cropped and downloaded data for 1 km 2 squares centered on the eddy-flux observation towers in both forests. The Copernicus Global Land Cover Layers (CGLS-LC100 collection 3; Buchhorn et al 2020) show that the 1 km 2 squares of both forests are homogeneously and thoroughly covered with closed, deciduous needleleaf forest, which is equivalent to larch forest.
The equation used to calculate the NDVI for each pixel in each 1 km 2 square was as follows: where R and NIR are surface reflectance from the red and near-infrared bands, respectively. The NDVI in pixels with low-quality data was determined and removed with bit values in QC_250m from MOD09GQ (bit0-1 was used for overall data quality) and state_1km from MOD09GA (bit0-1, bit2, and bit8-9 were used for cloud density, cloud shadow, and cirrus presence, respectively). We resampled MOD09GA data to a resolution of 250 m from 1 km using the nearest-neighbor method. Daily NDVI data were composited monthly, and monthly medians were calculated for each pixel. The monthly medians for 16 pixels per site were averaged spatially to produce a monthly mean value for each site. Thus, 100 monthly mean NDVI values were obtained for each forest (5 months yr −1 for 20 years). NDVI data for 2000-2019 for both forests are presented in figure 2(a).

Meteorological data
Because the tower-based in situ meteorological observations did not cover the entire study period, we obtained additional meteorological data captured

Data analysis
Mann-Kendall trend analysis, path analysis, and correlation analysis were used to examine monthly anomalies in NDVI, temperature, and precipitation over 20 years for both forests. Anomalies in May, June-August, and September were analyzed separately. In addition to the total precipitation for the month or season of interest, the total precipitation from May of the previous year to April of the current year (hereafter, previous May-April) was also included in the analyses because it is an important measure of soil water content in Siberian larch forests (Ohta et al 2008). Monthly averages for 2000-2019 were taken as climatological values in calculating anomalies. Mann-Kendall trend analyses were used to examine 20 year trends for the anomalies. Path analyses were performed to examine causal relationships between the NDVI, temperature, and precipitation. In this study, we focused on differences in the NDVI-temperature-precipitation relationship between the two forests. To this end, the chi-square statistic was used to determine the statistical significance of each model. Correlation analyses were used to evaluate the similarity of time courses of anomalies between the two forests.
In addition, we analyzed MCD43A4-derived NDVI data following the aforementioned protocol to evaluate the robustness of our analyses across different MODIS data products, because NDVI time series can be inconsistent, even between MODIS Terra and Aqua (Fensholt and Proud 2012), and satellite observations can be affected by the geometry of the sensor, the solar illumination angle, and the Earth's surface (Guay et al 2014). The MCD43A4 product includes daily reflectance derived from 16 days of MODIS Terra and Aqua data adjusted by a nadir bidirectional reflectance distribution function.
The analyses were conducted with CRAN R (R Core Team 2017) and associated R packages such as raster (Hijimans 2020), EnvStats (Millard 2013), stats (R Core Team 2017), and sem (Fox 2006). In the present study, p < 0.05 was taken to indicate statistical significance.

Time courses of anomalies
The time course of the NDVI anomalies over the past 20 years differed between Spasskaya Pad and Elgeeii (figure 3). At Spasskaya Pad, the most negative NDVI anomaly found for June-August was in 2003 (−0.11), whereas the most positive anomaly found was in 2001 (0.06). However, there was no apparent increasing or decreasing trend overall. At Elgeeii, the most negative NDVI anomaly found for June-August was also in 2003 (−0.14), whereas the most positive anomaly found was in 2017 (0.10), with the NDVI exhibiting a linearly increasing trend. NDVI anomalies in May ranged from −0.04 to 0.04 at Spasskaya Pad and from −0.08 to 0.13 at Elgeeii. In September, NDVI anomalies ranged from −0.09 to 0.19 at Spasskaya Pad and from −0.15 to 0.17 at Elgeeii. In addition, NDVI anomalies were not significantly correlated between Spasskaya Pad and Elgeeii for any season (p > 0.05). More apparent differences in the NDVI time series between the two forests were observed from analyzing NDVI trends over the past 20 years (table 1). Specifically, the NDVI exhibited a significant increasing trend for May and June-August at Elgeeii (p < 0.05), whereas no significant trends were observed for any season at Spasskaya Pad (p > 0.05).
The time series of temperature anomalies of the two forests were similar (figure 4), whereas the precipitation time series differed (figure 5). For both forests, the most negative temperature anomaly recorded for June-August was in 2004 (−6.6 • C and −5.9 • C for Spasskaya Pad and Elgeeii, respectively), whereas the most positive anomaly recorded was in 2008 (3.7 • C and 5.3 • C, respectively). Accordingly, temperature was strongly and positively correlated between the two forests (R 2 > 0.74, p < 0.01). Conversely, precipitation was not significantly correlated between the two forests (p > 0.05). At Spasskaya Pad, the most negative and most positive precipitation anomalies recorded for June-August were in 2001 (−66 mm) and 2006 (82 mm), respectively, whereas the respective years for Elgeeii were 2003 (−86 mm) and 2008 (57 mm). In contrast to the significant increasing NDVI trend, temperature and precipitation did not exhibit significant trends over the past 20 years for either forest (table 1; p > 0.05).

NDVI-meteorology relationship
Path analyses showed that the NDVI anomalies for May and June-August at Elgeeii were significantly associated with temperature and previous precipitation (p < 0.05), respectively, whereas no relationships were found between NDVI and temperature or precipitation anomalies at Spasskaya Pad (p > 0.05; figure 6). At Elgeeii, the NDVI in May was positively correlated with temperature (p < 0.01), whereas the NDVI in June-August was positively correlated with precipitation during the previous May-April period (p < 0.01). The path coefficient for the association between precipitation in June-August and NDVI in the same season was greater than that for temperature but not significant (p > 0.05). A significant chisquare test result was obtained only for the model that included NDVI anomalies in June-August at Elgeeii (p < 0.05), which implies that additional unknown factors play a role in NDVI anomalies in other months and in the Spasskaya Pad forest. Individual comparisons of NDVI anomalies and meteorological conditions are presented in figure S1 (available online at stacks.iop.org/ERL/17/025003/mmedia).
Detrended NDVI anomalies for May and June-August at Elgeeii were also associated (p < 0.05) with temperature and previous precipitation, respectively (figure 7). For June-August, precipitation in the previous May-April period was the only factor significantly associated with detrended NDVI anomalies, which implies that it was the primary controlling factor at Elgeeii. A significant chi-square result was obtained for the model that included detrended NDVI anomalies in June-August (p < 0.05). Thus, at least for June-August, each boreal forest had a different NDVI-meteorology relationship. However, the 20 year NDVI trend was likely controlled by factors that were not examined in this study, as no significant trend over the past 20 years was observed for temperature or precipitation (table 1).

Robustness of NDVI analyses based on MOD09GQ and MCD43A4 data
The NDVI data derived from the MODIS Terra MOD09GQ product were strongly correlated (R 2 = 0.875, p < 0.01) with those from MCD43A4, in particular for June-August ( figure 8(a)). The correlation for September was also significant (R 2 = 0.364, p < 0.01), whereas the one for May was not significant. Thus, MOD09GQ-derived NDVI data were similar to those from MCD43A4, especially for June-August, such that consistent results from the comparison of NDVI-meteorology relationships between the two forests were obtained ( figure 8(b)). For the Elgeeii forest, path coefficients indicating associations between temporal variation in NDVI anomalies and precipitation in the previous May-April period and June-August were significant.

Discussion
The difference between the two sites in the significance of NDVI trends over the period 2000-2019 (figure 3 and table 1) supports our first hypothesis. The observed annual rate of increase in NDVI for June-August at Elgeeii (0.004 yr −1 ) is comparable to the typical rate for terrestrial ecosystems in the pan-Arctic region (  Our second hypothesis is validated by the different NDVI−meteorology relationships between Spasskaya Pad and Elgeeii (figure 6). However, this difference was insufficient to explain the difference in NDVI trends between the two forests, as precipitation and temperature did not exhibit significant trends over the past 20 years (table 1). Moreover, significant relationships between the NDVI and meteorological conditions were found for the Elgeeii forest when detrended NDVI anomalies controlled by temperature in May and precipitation in the previous May-April period were used (figure 7). With MCD43A4derived NDVI data, precipitation during the present season was also a significant factor (figure 8). Therefore, in this study, both hypotheses stating that 20 year NDVI trends and NDVI-temperature-precipitation relationships would differ between larch forests with and without waterlogging-induced disturbances were validated. Nonetheless, further studies are needed to reliably determine linkages between spatially variable NDVI trends and vegetation-environment interactions in the Siberian boreal region.
A possible explanation for the less-apparent relationship between the NDVI and meteorological conditions at Spasskaya Pad (figure 6) is the year-to-year change in contribution from understory vegetation and canopy-forming trees after larch diebacks from waterlogging . These diebacks have resulted in understory vegetation contributing substantially to the overall features of Siberian larch forest ecosystems (Iida et al 2009, Suzuki et al 2011, Loranty et al 2018, Nagai et al 2020, with proportional contributions to total NDVI of 50% or more recorded in some cases (Suzuki et al 2011, Loranty et al 2018. Furthermore, the environmental responses of understory vegetation in boreal regions sometimes differ from those of canopy-forming trees (Ikawa et al 2015). In the Spasskaya Pad forest, understory vegetation was less vulnerable to excess wetting of the surface soil layers than were canopy-forming trees, according to carbon flux data . Thus, waterlogginginduced changes in forest composition, in particular the cover of understory vegetation and canopyforming trees, would likely affect the overall response of this forest ecosystem to environmental changes. The less apparent relationship between the NDVI and meteorological conditions at Spasskaya Pad may also be due to differences in geological origin giving rise to distinct forest ecological features at Spasskaya Pad and Elgeeii. Differences in interannual variation in precipitation were observed between Spasskaya Pad and Elgeeii (figure 5(c)), in addition to fundamental differences in vegetation biomass, precipitation amounts, and soil properties , Hiyama et al 2021. Thus, it is challenging to elucidate the factors underpinning the weak relationship between the NDVI and meteorological conditions at Spasskaya Pad.
The significant relationship between the NDVI and temperature or precipitation depending on the season at Elgeeii (figures 6-8) supports the wellknown assertion that both temperature and summer precipitation are important drivers of vegetation and carbon dynamics in Siberian boreal-forest ecosystems (Ohta et al 2008, Hiyama et al 2021. The positive correlation between the NDVI and temperature anomalies in May was likely due to warming-induced relief of the cold environment, as temperature fundamentally regulates the growth of vegetation in the pan-Arctic region (Chapin et al 2010, Zhu et al 2016. The importance of precipitation, in particular that in the previous May-April period, for the NDVI in June-August agrees with the result of previous tree ring analyses indicating a significant dependence of Siberian larch tree growth on precipitation in the previous year (Rigling et al 2001, Tei et al 2014, Tei and Sugimoto 2020. In addition, Ohta et al (2008) reported apparent relationships between summer precipitation in the previous year and soil water content in the current year. Hence, although NDVI trends were not strongly associated with meteorological conditions (figures 6, 7 and table 1), the overall change in summer precipitation is increasingly recognized as an important driver of vegetation and carbon dynamics in Siberian boreal-forest ecosystems (Hiyama et al 2016). Figure 6. Path models describing causal relationships between the NDVI and meteorological metrics such as air temperature and current-and previous-year precipitation, including possibly covarying relationships among these three meteorological factors. Statistically significant path coefficients are indicated with asterisks ( * and * * , p < 0.05 and p < 0.01, respectively). The anomalies presented in figures 3-5 were used in these path models.
One possible driver of the increasing NDVI trend at Elgeeii is the increasing concentration of atmospheric CO 2 , namely, CO 2 fertilization (Norby et al 2005, Ueyama et al 2016, Zhu et al 2016. On a global scale, CO 2 fertilization at a rate of 2.2 ppm yr −1 in recent decades has likely increased gross primary production by 1.2%, equivalent to 60% of the total increase in global gross primary production from 2000 to 2014 (Ueyama et al 2020). Similarly, Zhu et al (2016) estimated that CO 2 fertilization has contributed to 70% of the total increase in the LAI on a global scale. In boreal forest ecosystems, CO 2 fertilization was estimated to stimulate gross primary production by 2%-4% during the 2002-2014 period (Ueyama et al 2016), although the warming climate is the most important driver in this region (Zhu et al 2016). Thus, the increasing CO 2 concentration is likely the most important factor among those unexamined in the present study driving an increase in the NDVI at Elgeeii. Path model describing causal relationships between detrended NDVI anomalies and meteorological metrics such as air temperature and current-and previous-year precipitation, including possibly covarying relationships among these three meteorological factors, in May and June-August (JJA) for the Elgeeii forest. Statistically significant path coefficients are indicated with asterisks ( * and * * , p < 0.05 and p < 0.01, respectively). The anomalies presented in figures 3-5 were used in this path model.

Conclusions
We found apparent differences in NDVI trends for the 2000-2019 period between two boreal larch forests with and without waterlogging-induced disturbances between 2005 and 2008. Relationships between the NDVI and meteorological anomalies also differed between the two forests, although this difference was insufficient to explain the dissimilarity in NDVI trends. Based on the findings, we infer that Siberian boreal larch forests with nonsignificant NDVI trends might have experienced waterlogging-induced dieback of larch trees and a concomitant increase in water-resistant understory vegetation. Future studies should consider the possibility of ecological disturbance and/or changes in the composition of vegetation, even in boreal-forest ecosystems that exhibit no significant NDVI trends in recent decades.