Different Temporal Stability and Responses to Droughts between Needleleaf Forests and Broadleaf Forests in North China during 2001–2018

Abstract

Droughts can affect the physiological activity of trees, damage tissues, and even trigger mortality, yet the response of different forest types to drought at the decadal time scale remains uncertain. In this study, we used two remote sensing-based vegetation products, the MODIS enhanced vegetation index (EVI) and MODIS gross primary productivity (GPP), to explore the temporal stability of deciduous needleleaf forests (DNFs) and deciduous broadleaf forests (DBFs) in droughts and their legacy effects in North China from 2001 to 2018. The results of both products showed that the temporal stability of DBFs was consistently much higher than that of DNFs, even though the DBFs experienced extreme droughts and the DNFs did not. The DBFs also exhibited similar patterns in their legacy effects from droughts, with these effects extending up to 4 years after the droughts. These results indicate that DBFs have been better acclimated to drought events in North China. Furthermore, the results suggest that the GPP was more sensitive to water variability than EVI. These findings will be helpful for forest modeling, management, and conservation.

1. Introduction

Forests are important ecosystems globally, providing nearly half of the total gross primary production for terrestrial ecosystems [1]. Yet, there is growing evidence that a drought could trigger [2,3,4] or even exacerbate [5] the mortality of trees, thus weakening the ecosystem functions of forests [6]. At the global scale, there has been an increasing trend in the frequency of extreme drought events in recent years [7]. However, due to spatial heterogeneity, differences in the spatial and temporal characteristics of droughts at the regional scale [8] may cause variable responses from forests to these extreme events.

How forests might respond to droughts is a critical issue with many uncertainties. From an ecological perspective, the response of forest ecosystems to droughts can be characterized by many aspects, such as temporal stability, resistance, and resilience. These aspects are all part of the overall ecosystem stability [9]. Research has shown that the ecosystem stability against droughts may vary among different regions [10], forest types [11], and even tree species [12].

Forest responses to droughts may also be highly influenced through legacy effects. Numerous studies have shown that tree growth might not respond to a drought immediately, as there is a lag period between a drought and the forest response that can range from several months to as long as several years [13,14,15]. Studies have suggested differences in this time lag based on forest types and tree species, but a quantitative comparison is still needed.

Using tree-ring data is one of the most popular methods for studying forest stability [16,17]. As a field measurement, tree-ring records can show from decades to hundreds of years of tree growth. These long records reflect multiple environmental conditions over many years. Thus, scientists have used tree-ring data to reconstruct historical climates [18]. However, in situ measurements such as tree-ring records are difficult at a regional scale. Even the world’s largest public tree-ring database, the International Tree-Ring Data Bank (ITRDB) [14], lacks records in many regions, including North China. As a result, studies about forest stability against droughts at a global or hemispheric scale might lack data for forests in this region [2,16]. As such, remote sensing data sources serve as another major tool for studying forest stability. Although remote sensing datasets do not cover as long a period as the tree-ring records, they are good compliments for decadal analyses, which can cover the whole globe [19]. Further, remote sensing-based vegetation indices have been used to quantify resilience across biomes [11], to measure leaf and canopy level activities [20], and to study the responses of forests to droughts at a regional scale [21]. The key limitation with remote sensing data is that they only record the spectral information; thus, they are generally an indirect measurement of the structural and physical characteristics of vegetation. Therefore, the vegetation characteristics reflected by different remote sensing products need to be carefully considered when using these datasets. There is still an outstanding need for research in this area to be fully developed.

North China has experienced several severe-to-extreme drought events in the last two decades [22,23]. To the best of our knowledge, research that has been conducted on the forests in this area includes using tree-ring data to study the effect of historical extreme droughts on forests [24], the response of the forest phenology to preseason droughts [25], and the influence of forest coverage during droughts [21]. However, ecosystem stability for different forest types and their response time lag to droughts remain unclear. Using two remote sensing-based datasets and one meteorological drought index, we tried to answer two questions: (1) was there any difference between the temporal stability of deciduous broadleaf forests and deciduous needleleaf forests in North China during 2001–2018, and (2) how long did it take the two types of forests to respond to droughts? The answers to these questions provide more information for quantifying ecosystem stability and clarifying the impact of the forest type in response to droughts. These results benefit forest modeling, management, and conservation.

2. Materials and Methods

2.1. Study Area

North China covers an area of about 1,515,650 km². The southeast region of North China is in the temperate monsoon ecozone, while the northwest part of this area is in the semi-arid ecozone. Despite the climate heterogeneity of this region, forests are distributed broadly in this area and provide important ecosystem services (Figure 1b, the colored pixels). Typical forest types in North China include deciduous broadleaf forests (DBFs), deciduous needleleaf forests (DNFs), evergreen needleleaf forests (ENFs), and mixed forests. Over the past two decades, forest cover has significantly changed in this area, mostly because of deforestation, afforestation, and other human activities. This land cover change might lead to a more complex spatial pattern of the response of forests to droughts. To disentangle the responses of various forests to droughts without land cover change, we focused on the forests which were constant during 2001–2018 in North China. The typical forest types constant during this period were DBFs and DNFs (see Section 2.2.4 for detailed selection criteria and procedures).

2.2. Data

2.2.1. Remote Sensing-Based Vegetation Growth Products

Two remote sensing-based vegetation growth products were used in our study: the enhanced vegetation index (EVI) and gross primary productivity (GPP). The EVI has been widely used in ecological research [26,27,28], and furthermore, compared with the commonly used normalized difference vegetation index (NDVI), it has improved performance over high biomass regions such as forests [29] and reduces the influence of water vapor and clouds. We used the MODIS 500 m 16-day EVI (MOD13A1 Version 6, https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD13A1, accessed on 26 March 2021) here. A MODIS GPP product (MOD17A2H collection6, https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD17A2H, accessed on 27 March 2021) with a 500-m spatial resolution and 8-day temporal resolution was employed as the GPP data source. This data was generated from an improved light use efficiency (LUE) model with the MODIS FPAR and PAR data [30].

Both of the MODIS datasets were mosaiced and extracted via the Google Earth Engine (GEE) platform [31]. To make the time resolution comparable to the drought index, we recalculated the datasets for each month using the days as weights. For example, the first record of the EVI in 2001 was on 1 January; the second record was 16 days later, on 17 January; and the third record was on 2 February. To obtain the monthly record, we defined the day weight of the first record as 16/31 in January, the second record as 15/31 in January, and the third as 1/28 in February in a non-leap year (1/29 in February in a leap year). This is also an official processing procedure in the MODIS group for producing the monthly NDVI [32].

2.2.2. Drought Index

To detect drought events in our study area, we used the Standardized Precipitation-Evapotranspiration Index (SPEI). The SPEI is a state-of-the-art meteorological drought index. It considers not only the temperature and precipitation (PPT) as the traditional drought index does (e.g., PDSI and SPI) but also potential evapotranspiration (PET), which accounts for the real water balance for terrestrial ecosystems [33]. In this study, we downloaded the 12-month windowed SPEI dataset from the global SPEI database (https://spei.csic.es/database.html, accessed on 18 January 2021), since SPEI-12 is a good representation for annual water balance [17,34] and forests are sensitive to annual drought events [35,36]. The spatial resolution of this dataset was 0.5 degrees. To align the pixels with the MODIS EVI and GPP for spatial statistical analyses, we resampled the SPEI pixel into 500 m following the nearest neighbor rule.

Generally, the levels of the water balance conditions could be classified by the thresholds of the SPEI as follows (

Table 1. Thresholds of Standardized Precipitation-Evapotranspiration Index (SPEI) values for water conditions.

Standardized Precipitation-Evapotranspiration Index (SPEI) Value	Condition
>2	Extremely wet
(1.5, 2]	Severely wet
(1, 1.5]	Moderately wet
(0, 1]	Normally wet
[−1, 0)	Normal drought
[−1.5, −1)	Moderate drought
[−2, −1.5)	Severe drought
<−2	Extreme drought

).

2.2.3. Climate Data

Climate factors significantly influence vegetation activities. To analyze the correlation between vegetation growth (EVI and GPP) and the drought condition (SPEI), we also included other potential dominant climate factors, specifically long-term temperature (TMP) and shortwave radiation (RAD), in our study. Since the SPEI dataset used the CRU climate dataset as an input, we chose the CRU 0.5 degrees monthly dataset (version TS 4.03, https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.03/, accessed on 21 September 2019) as the temperature data source. Radiation data were not included in the CRU database, so instead, we downloaded the monthly downward solar radiation from the state-of-the-art NCEP/NCAR Reanalysis II dataset (ftp://ftp2.psl.noaa.gov/Projects/Datasets/ncep.reanalysis2/gaussian_grid, accessed on 26 July 2021). These RAD data were originally binned into the Global T62 Gaussian grid, which equated to approximately 1.85 longitude degrees, and it covered the period from 1 January 1979 to 31 July 2021. Similar to the SPEI data, we selected the CRU-TMP data and the NCEP/NCAR-RAD data during 2001–2018 and resampled the data into a 500-m spatial resolution for our analysis.

2.2.4. Vegetation Category

The Vegetation Map of China (vector format, 1:1,000,000) from the Institute of Botany of the Chinese Academy of Sciences [37], as well as the MODIS landcover product MCD12Q1, were used for extracting and categorizing different forest types. The Vegetation Map of China is derived from multi-source information, including forest inventory data, plot-level measurements, and remote sensing data. The detailed base data help to bring the Vegetation Map of China closer to the actual vegetation distribution. Unfortunately, this map is generated with data over many years and thus may not reflect land cover change. Meanwhile, the MCD12Q1 is a dynamic dataset recording the landcover type for each year from 2001 to 2019 at a 500-m spatial resolution. Thus, we converted the Vegetation Map of China into 500-m resolution raster data and compared them with the MCD12Q1 for each year and pixel. We then chose the pixels which were constant for all years in the MCD12Q1 and consistent between the Vegetation Map of China and the MCD12Q1 for our study. As a result, most of the constant forest cover from 2001 to 2018 was deciduous broadleaf forests (DBFs, 50,737 pixels) or deciduous needleleaf forests (DNFs, 26,425 pixels). Only 822 pixels were categorized as constant evergreen needleleaf forests, and neither constant nor inconstant evergreen broadleaf forest pixels were found in the study area. Hence, our research would focus on the DNF and DBF biomes in North China.

2.2.5. Temporal Stability

To quantify the temporal stability of the forests in North China during 2001–2018, we used a common definition [9,38] as follows:

$$\mathrm{Temporal}\mathrm{Stability}\left(\mathrm{TS}\right)=\frac{ \mu }{\sigma }=\frac{\overline{GS}}{\sqrt{\frac{1}{n-1}{\displaystyle {{\displaystyle \sum }}_{i=1}^n}\left(G{S}_i-\overline{GS}\right)}}$$

(1)

where

$$\overline{GS}=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}G{S}_i$$

(2)

for the EVI,

$$G{S}_i=\frac{1}{5}{{\displaystyle \sum }}_{j=5}^{j=9}EV{I}_j$$

(3)

for the GPP,

$$G{S}_i={{\displaystyle \sum }}_{j=5}^{j=9}GP{P}_j$$

(4)

where μ and σ are the mean and the standard deviation of the growing season vegetation growth condition $GS$ across year i (i = 2001, 2002, … 2018), respectively, and j stands for the jth month within a certain year. For the EVI, the $G{S}_i$ is the average of monthly EVI; for GPP, the $G{S}_i$ is the sum of monthly GPP of the growing season.

We defined the growing season for both DNFs and DBFs as being from May to September because previous remote sensing studies have shown that the growing season of both DNFs and DBFs in North China covered from May to September for the past two decades [39].

For the spatiotemporal analysis, we retrieved the values of the EVI and GPP for each growing season from 2001 to 2018 to calculate the TS_EVI and TS_GPP for each pixel.

2.3. Statistical Methods

2.3.1. Interannual Trend Detection

To detect the temporal dynamics of forest growth and water balance in our study area from 2001 to 2018, we deseasonalized the monthly EVI, GPP, and SPEI and then used Ordinary Least Squares (OLS) regression to calculate the interannual trend.

For time series analysis, the records across continuous time generally contained characteristics about a seasonal cycle, interannual change (trend), variation, and noise. It was necessary to remove the seasonal component before calculating the trend. Here, the seasonal cycle was defined as the average value of multiyear records in each month from January to December. Deseasonalization would remove this cycle and leave only the interannual trend and variation. The “deseason” function in the climate data toolbox (CDT) for Matlab was used for the deseasonalization. After that, the Matlab function “regress” was used to analyze the interannual trend of the EVI, GPP, and SPEI based on the Ordinary Least Squares (OLS) regression method. The deseasonalized time series were only used for describing the interannual trends of EVI, GPP, and SPEI in this study.

2.3.2. Partial Correlation Analysis

Climate factors such as temperature, water availability, and radiation were considered to be the most important driving factors of ecosystem stability. To quantify the correlation between the temporal stability and the water balance while eliminating the colinear influence of other climate factors, we utilized partial correlation analysis. We set the mean value of the growing season TMP and RAD as controlling variables and then calculated the partial correlation coefficient between TS and the mean value of the growing season SPEI for all valid pixels. This analysis was performed with the Matlab function “partialcorr”. As a supplement, we also calculated the partial correlation coefficients of the growing season between TMP and TS with controlled RAD and SPEI and the partial correlation coefficients of the growing season between RAD and TS with controlled TMP and SPEI, respectively. The complete results are attached to Appendix A as

Table A1. Partial correlation coefficients of SPEI_SD, temperature and radiation with TS_EVI and TS_GPP. In this study, the p values of all results were less than 0.001.

	TSEVI			TSGPP
	All	DNF	DBF	All	DNF	DBF
SPEISD	−0.24	−0.20	0.01	−0.05	−0.23	−0.09
TMP	−0.29	0.06	−0.42	−0.27	0.03	−0.35
RAD	0.05	−0.25	0.22	0.07	0.19	0.16

.

2.3.3. Cross-Correlation Analysis

Cross-correlation analysis is an efficient way to detect a time lag between two signals. In this study, we used this method to detect if forest growth during the growing season had a legacy effect on the SPEI. For example, we first ran Pearson correlation analysis between the EVI during May–September from 2001 to 2018 and the SPEI during the same period, which implied a 0-month lag, or no legacy effect. Secondly, we moved the SPEI series one month earlier, or April–August from 2001 to 2018, and analyzed its correlation with the EVI during May–September from 2001 to 2018, which implied a 1-month lag of the forest growth (Figure 2). Since studies have demonstrated that global forests generally have 1–4 years of lag to droughts, here, we set the shifted lag step from 0 months to 60 months (i.e., 5 years) and calculated the correlation coefficient (r) for each valid pixel. After that, we detected the r value of the largest absolute value. The corresponding lag step corresponded to the length of the legacy effect for each forest pixel.

3. Results

3.1. The Interannual Variation of Water Balance and Forest Growth from 2001 to 2018

We extracted the monthly SPEI12, EVI and GPP values for each constant forest pixel and then looked at the spatial average to draw the time series. As Figure 3a shows, the water balance for forests in North China significantly increased from 2001 to 2018 (slope = 0.003, p = 0.0011), which suggested an increase in water availability for this region. However, despite this trend, forests in this area were generally under water stress, since SPEI12 was negative for most months. The mean SPEI12 for 2001–2018 was −0.7191, which suggested slight droughts. Furthermore, during these 18 years, there were four severe drought events (2003: SPEI_min = −1.59; 2007: SPEI_min = −1.81; 2008: SPEI_min = −1.98; and 2017–2018: SPEI_min = −1.57).

The EVI (slope = 0.0001, p = 0, Figure 3b) and GPP (slope = 0.015 g C/m², p = 0, Figure 3c) increased significantly as well over this time period. These results suggest that increased water availability influenced forest growth, but this impact might not have been linear or immediate. Moreover, the fluctuations in the summer values of the EVI and GPP in the four severe drought years were different. Although both the deseasonalized EVI and GPP had much lower values in summer in 2003 compared with the general season cycle, the EVI was higher than normal in 2007, an extreme drought year. Meanwhile, the deseasonalized GPP did not increase in the summer of 2007. Similar patterns were observed in 2017 (Figure 3b,c). These discrepancies might be explained by different physical characteristics reflected by the two indices, which influenced their response times to droughts.

3.2. The Temporal Stability of Forest Growth in North China from 2001 to 2018

Though the TS_EVI and TS_GPP showed similar distributions, the magnitudes of the temporal stabilities generated from these two datasets were different (Figure 4a,c versus Figure 4b,d). The range of TS_EVI was [6.43, 45.81], and the range of TS_GPP was [0.96, 27.94]. The higher range of TS_EVI implied that the temporal variations of the EVI and GPP were out of sync. Because the constant forest pixels were distributed sparsely in the study area and without an apparent pattern, the spatial distribution of the TS are shown in the appendix (Figure A1). To determine which component part (mean or standard deviation, SD) was the dominant factor of the TS, the mean and SD of the EVI and GPP were plotted against the TS (Figure 4). We found that the standard deviation, which stands for the variation of interannual change, explained more of the TS values.

3.3. Partial Correlation between TS and SPEI

Since the TS_EVI and TS_GPP were both highly correlated with the SD of the mean GS SPEI, we analyzed the partial correlation between the TS_EVI and TS_GPP and the SD of the mean GS SPEI with the mean growing season temperature (TMP) and mean growing season shortwave radiation (RAD) controlled. The correlation between the SPEI and TS_EVI was similar to that between the SPEI and TS_GPP among the DNFs with a negative result. However, the EVI and GPP for DBFs had positive results as well as a weak correlation with the water balance (Figure 5). The r value from both the EVI and GPP for DBFs was less than 0.1. This could explain the significant difference between TS_EVI and TS_GPP when the TS values of all the valid pixels were analyzed together. These results suggest that the canopy greenness and photosynthetic rates responded to droughts differently. Furthermore, DNFs were more sensitive to changes in the water balance than DBFs. We also calculated the partial correlation coefficients of TMP and RAD for reference. The complete results are shown in Table A1.

3.4. Legacy Effect of the Water Condition on Forest Growth

The sensitivity of forest growth to the water balance can influence an ecosystem’s stability [15]. Additionally, as numerous studies have shown, legacy effects can also influence the response of vegetation growth to shifts in the water balance [20,40,41], as well as fluctuations in the time series of the EVI and GPP. Since the EVI and GPP were asynchronous in this study (Figure 3), a further exploration with cross-correlation analysis was conducted between the two vegetation growth condition data and SPEI12. Considering research demonstrating the 1–4 years’ legacy effect of forests with droughts [13,15], we set the maximum lag threshold as 60 months (5 years).

As a result, the SPEI showed key legacy effects on both the growing season EVI and GPP of forests pervasively (Figure 6). Although the legacy effects varied from 0 to 60 months, the pixel number binned to each lag length followed a right-skewed distribution for the EVI (Figure 6a) and a left-skewed distribution for the GPP (Figure 6b). For the growing season EVI, 85.7% of the forest pixels showed a lag within 48 months. For the growing season GPP, much shorter lag responses were exhibited. However, 26.9% of the forest pixels showed only a 1-month lag to the SPEI (20,734 out of 77,162 pixels). Most of the correlations were positive, with a proportion of 81.6% from the EVI and 77.2% from the GPP.

When the legacy effect results were categorized with DNFs and DBFs, the DBFs’ response characteristics were more consistent than the DNFs. For DBFs, most of the pixels illustrated a lag from 3 to 4 years for the EVI to the SPEI and a 1-month lag for the GPP to the SPEI. However, the histograms of the DNFs were more alike in the multimodal patterns (Figure 6).

When significance was considered, the patterns were similar (Figure 6c,d). This huge discrepancy implied that the EVI generally responded to shifts in water availability much faster than the GPP, especially for the DBFs.

4. Discussion

4.1. Stability Comparison between DBFs and DNFs

We found no significant difference between the temporal stability of DBFs and DNFs (Figure 4); the range and distribution were similar. However, when the correlations among the SPEI, TMP, and RAD with TS were disentangled, the results suggested that the TS of the two forest types responded to the water balance and other climate factors differently (Figure 5 and Table A1). For DNFs, water balance was an important factor. The smaller the variation of the interannual SPEI was, the more stable the DNFs would be (Table A1). In contrast, for DBFs, temperature was the limiting factor. High temperatures stressed DBFs’ growth. However, shortwave radiation was similarly important for both DNFs and DBFs. Hence, it is possible that ecological factors connected with latitude (e.g., solar radiation and temperature) are impact factors for the temporal stability of both forest biomes, and water may only influence certain forest types. Research based on global remote sensing data and random forest model simulations support this hypothesis. Previous results have illustrated that the mean annual shortwave radiation, mean annual temperature, and precipitation contribute the most to ecosystem stability [42].

Interestingly, even though the region between 40° N and 50° N of North China experienced extreme droughts, the DBFs between the two latitudes still maintained a relatively high level of temporal stability (Figure A1 and Figure A2). This could be explained by the optimal temperature threshold, because the optimal photosynthesis temperature for DBFs was generally higher than the observed max air temperature [43,44], and droughts caused by high temperatures in this area did not significantly decrease the photosynthesis rate.

Furthermore, the EVI is a spectral index that reflects the pigment and cell structure characteristics of leaves at the canopy scale. The GPP is the gross photosynthesis rate of an ecosystem. The two datasets stand for different perspectives of vegetation growth, and thus the essential differences between optical and physiological characteristics could be used to explain the discrepancy between the temporal stability from the EVI and GPP. Previous research has shown that the canopy structure of a broadleaf forest was not as sensitive to drought as the xylem formation [42,45]. Experimental manipulations also have demonstrated that several forest species could acclimate to warming by reducing their respiration at a year-to-decades time scale [46], which partly explains the unbelievable stability of DBFs during extreme droughts and the weak positive correlation between the temporal stability and SD of the SPEI. Our future work will further explore the contributions of temperature, radiation, and water to ecosystem stability through droughts.

4.2. Legacy Effect Comparison between DBFs and DNFs

The DBFs not only showed higher temporal stability than DNFs, but they also demonstrated more similarity in the lengths of their legacy effects with droughts. This result was consistent with a previous study which found a longer legacy effect for DNFs than DBFs [19]. Our results based on the EVI and GPP might also be a bridge for previous research about drought legacy effects on forests based on different data sources. Some studies showed a legacy length less than 4 years [15,19], while others showed a legacy no longer than 2 years [13]. The former research works were based on tree-ring data across Europe and America, which reflected the xylem accumulation or carbon partitioning strategy. The latter research was based on multi-source GPP data, which reflected ecosystem photosynthesis or carbon fixation. The diverse time lags implied that forests respond to water differently under various ecological processes. These differences call for more field measurements and more research about carbon partitioning strategies during droughts. Moreover, the relatively converging long time lag of DBFs suggests that different broadleaf forest communities have acclimated to droughts. This result provides a reference in studying the response of photosynthesis to droughts in different vegetation types.

5. Conclusions

We explored the temporal stability of DBFs and DNFs in North China with two remote sensing-based datasets, MODIS EVI and GPP, by quantifying their correlation with climate factors and legacy effects. Our results suggest that the two datasets were able to depict similar spatial distributions of temporal stability for forests in North China. These datasets were also able to capture the long-lasting legacy effects of forest responses to the SPEI. The results suggested that the remote sensing-based EVI and GPP could be used as substitutes for detecting responses of forests to droughts in areas lacking field measurement records. Meanwhile, the difference among the legacy effect detected from different datasets reminds us to select data according to the purpose of the analysis. Despite the multiple data sources, divergent temporal stabilities and legacy effects between DNFs and DBFs suggested that DBFs were more acclimated to droughts. Our further work will consider disentangling the contribution of climate factors to forest stability.