The Impact of Seasonality and Response Period on Qualifying the Relationship between Ecosystem Productivity and Climatic Factors over the Eurasian Steppe

Abstract

As climate change intensifies, surface vegetation productivity and carbon exchange between terrestrial ecosystems and the atmosphere are significantly affected by the variation of climatic factors. Due to the sensitivity of grasslands to these climatic factors, it is crucial to understand the response of vegetation greenness, or carbon exchange within grasslands, to environment factor dynamics. In this study, we used solar-induced chlorophyll fluorescence (SIF), precipitation (P), vapor pressure deficit (VPD), evaporative stress (ES), and root zone soil moisture (RSM) derived from remote sensing, reanalysis, and assimilation datasets to explore the response of vegetation greenness within Eurasian Steppe to climatic factors. Our results indicated deseasonlization based on the Seasonal-Trend decomposition using Loess (STL) method, which was an effective means to remove the seasonality disturbances that affect the qualification of the relationship between SIF and the four climatic factors. The response of SIF had a time lag effect on these climatic factors, and the longer the response period, the greater the impact on the correlation of SIF with P, VPD, ES, and RSM. We also found, among the four factors, that the response of SIF to ES was the timeliest. The findings of this study emphasized the impact of the seasonality and time lag effect on the dynamic response between variables, and provided references to the attribution and monitoring of vegetation greenness and ecosystem productivity.

1. Introduction

Irregular changes of environmental factors caused by climate change have significant impacts on the change of surface vegetation growth, resulting in erratic carbon exchange between terrestrial ecosystem and the atmosphere [1]. Grasslands are a potentially important sink for carbon in the terrestrial biosphere because of the high levels of carbon accrual and sequestration below ground in the soil, and up to 90% of grassland biomass is below ground. The levels of soil carbon are higher in grasslands than in forests, agroecosystems, or other ecosystems [2,3]. Carbon cycling within grasslands is sensitive to climatic factors [4] and vulnerable to climate change [5,6,7]. As climate change intensifies, it is crucial to understand the response of vegetation greenness or carbon exchange within grassland to environment factors’ dynamics.

Normally, the carbon exchange between terrestrial ecosystem and the atmosphere is quantified by gross primary productivity (GPP) and ecosystem respiration (RE) [8,9]. The solar-induced chlorophyll fluorescence (SIF), which is an optical signal emitted in the spectral range 650–850 nm from the chlorophyll a molecule in vegetation [10], demonstrates a tight relationship with photosynthesis at the ecosystem scale through ground-based and airborne hyperspectral measurements [11,12,13]. The global SIF observations from satellites have become available in recent years, which has led to tremendous developments in monitoring terrestrial photosynthesis and ecosystem productivity (i.e., GPP) at various spatial and temporal scales, which caused SIF to be increasingly used as a proxy for GPP on large scales [14,15,16]. Additionally, SIF was utilized to characterize vegetation greenness [17], determine the vegetation growing season [18], estimate evapotranspiration [19], and the like. Climatic factors, such as precipitation, temperature, solar radiation, transpiration, evaporation, soil moisture, etc., influence land-atmosphere carbon exchange through their impact on vegetation growth and greenness. Precipitation is the primary connection in the water cycle that provides for the delivery of atmospheric water to the Earth, and it controls ecosystem dynamics and drives biological processes [20,21]. The appropriate temperature and solar radiation, which are usually used as input variables in the calculation of some indices, are the elements that ensure the normal growth of vegetation [22]. The vapor pressure deficit (VPD) is calculated from the near-surface temperature and relative humidity, which can characterize the water stress in the air [23,24]. The soil moisture provides available water to vegetation roots and represents the water stress from the soil [25,26]. Evapotranspiration is the sum of transpiration from vegetation and the evaporation from the soil and land surface, and the important connection in the water cycle between the atmosphere and the land [27,28,29].

A clear understanding of the responses of vegetation productivity or ecosystem carbon uptake to the above climatic factors is helpful to extreme events monitoring, mechanism model establishment, etc. Previous studies evaluated the vegetation response to changes in air temperature [30,31], precipitation [32,33,34], and solar radiation [35,36], but fewer climatic variables were involved. However, since the response periods of vegetation to different climatic factors are diverse, it is also necessary to concurrently consider the relationships between grasslands and multiple climatic factors [37].

Here, we used SIF, precipitation (P), VPD, evaporative stress (ES), and root-zone soil moisture (RSM) derived from remote sensing, reanalysis, and assimilation datasets, to explore the responses of vegetation growth and productivity to climatic factors over the Eurasian Steppe. Before this, the deseasonlization method based on the Seasonal and Trend decomposition using Loess (STL) and Pearson’s correlation coefficient were utilized to probe the impact of seasonality/periodicity and the lag phases on qualifying the relationship between SIF and the four climatic factors. These two methods were also used to remove the disturbance of periodicity and to determine the response period between SIF and climatic factors.

2. Materials and Methods

2.1. Study Area

Grasslands, which are among the most prominent biogeographic regions on the Earth for human welfare and livestock, cover a wide range of vegetation types and geographic locations [7,38,39,40]. Eurasia contains nearly double the amount of total grassland ecoregions of the other continents [41]. The grasslands (27–125°E, 27–57°N, Figure 1) in Eurasia were thus selected as the study area, and it is the largest continuously continental grasslands in the world. The region has a complex topography such as lowlands, large plains, and high mountains, and is dominated by a temperate continental climate with hot summers and cold winters [42]. In addition, the precipitation in this area varies greatly from the northwest to the southeast. In recent decades, the grasslands are experiencing degradation and desertification caused by agriculture activities and disasters, which led to carbon loss from the grasslands. For degraded soils within the grassland, recovery is slow, taking 50 or more years [2]. These all result in damage to rangelands that are essential for animal husbandry development, regional ecological security, and sustainability [43].

2.2. Data

2.2.1. Land Cover

The Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Climate Modeling Grid Product (MCD12C1) was created using the supervised classification of MODIS reflectance data, and provides maps of the International Geosphere-Biosphere Program (IGBP), the University of Maryland (UMD), and the leaf area index (LAI) schemes at a 0.05° spatial resolution per year [44]. In this study, we followed the IGBP classification principle to extract the grasslands in the study area.

2.2.2. Solar-Induced Chlorophyll Fluorescence (SIF)

Solar-induced fluorescence (SIF) is mechanistically linked to photosynthesis and is shown to have a near-linear relationship with ecosystem GPP at the ecosystem scale [45,46]. In this study, the global ‘OCO-2’ SIF (GOSIF) dataset was used as the indicator of ecosystem productivity (i.e., GPP). The GOSIF was developed by Li and Xiao [13] using discrete OCO-2 SIF soundings, remote sensing data from the moderate resolution imaging spectroradiometer (MODIS), and meteorological reanalysis data through Cubist regression tree model. The GOSIF has a 0.05° spatial resolution and an 8-day and monthly temporal resolution. This product was available from March 2000 to December 2020.

2.2.3. Precipitation (P)

The Global Precipitation Measurement (GPM) precipitation dataset, which is organized as an international project led by the National Aeronautics and Space Administration (NASA) and the Japanese Space Agency (JAXA), is an international network of satellites that provide next-generation global observations of rain and snow [47]. Building upon the success of the Tropical Rainfall Measuring Mission (TRMM), the GPM concept centers on the deployment of a “Core” satellite carrying an advanced radar/radiometer system to measure precipitation from space and serve as a reference standard to unify precipitation measurements from a constellation of research and operational satellites. The Integrated Multi-satellitE Retrievals for GPM (IMERG) Level3 V06 monthly “Final” precipitation product with a 0.1° spatial resolution was used as a precipitation variable for climatic factors.

2.2.4. Air Temperature (T) and Relative Humidity (RH)

In this study, we used air temperature and relative humidity derived from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA 5) to calculate the VPD. The air temperature within the ERA 5 is the temperature of the air at 2 m above the surface of land, sea, or in-land waters, and calculated by interpolating between the lowest model level and the Earth’s surface, taking account of the atmospheric conditions. The relative humidity is the ratio of the partial pressure of water vapor to the equilibrium vapor pressure of water at the same temperature near the surface. Both have a 0.25° spatial resolution and monthly temporal resolution [48].

2.2.5. Evaporative Stress (ES) and Root Zone Soil Moisture (RSM)

The Global Land Evaporation Amsterdam Model (GLEAM) is designed to estimate the land surface evaporation and root zone soil moisture from remote sensing observations and reanalysis data [49,50]. Specifically, the Priestley-Taylor equation is used to calculate potential evaporation within the GLEAM based on the near-surface temperature and net radiation, while the root zone soil moisture is obtained from a multilayer water balance driven by precipitation observations and updated with the microwave soil moisture estimates [50]. The actual evaporation is estimated by constraining potential evaporation with a multiplicative evaporative stress factor based on the root zone soil moisture and vegetation optical depth estimates. The GLEAM version 3.5a provides the global daily and monthly potential and actual evaporation, evaporative stress conditions, and the root zone soil moisture from 1980 to 2020 at a spatial resolution of 0.25° (http://www.gleam.eu, accessed on 20 May 2021). In this study, the GLEAM v3.5a root zone soil moisture and evaporative stress were used.

The time period of SIF, P, T, RH, ES, and RSM we used in this study were from January 2001 to December 2020, and the 2001–2019 MCD12C1 land cover product was used. With regard to the spatial resolution, through the resample method which was ‘nearest’ in ArcGIS (https://desktop.arcgis.com/zh-cn/arcmap/10.6/map/main/mapping-and-visualization-in-arcgis-for-desktop.htm, accessed on 1 June 2021) version 10.6, we resampled the GPM precipitation (0.1°) and the GOSIF (0.05°) to 0.25°, that is, the same as the ERA5 and GLEAM datasets, which also meant the resolution of the pixel used in the next study and analysis was 0.25° × 0.25°.

2.3. Methods

2.3.1. Seasonal-Trend Decomposition Using Loess (STL)

The STL can decompose a time-series into additive components (trend, seasonality/periodicity, and the remainder) of variation by the application of loess smoothing models [51,52]. Further information on the method and parameters can be found in the original paper describing the STL method [53]. In this study, we used the STL to extract the seasonality of the time series of the variables mentioned (SIF, P, VPD, RSM, and ES) across each pixel over the study area. Then, we obtained the deseasonalized time series of this variables through the use of the observed time series minus the seasonal time series. The goal was achieved by Python 3.7 (https://www.python.org/, accessed on 4 June 2021) and the statsmodels library (https://www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.STL.html, accessed on 4 June 2021).

2.3.2. Pearson’ s Correlation Coefficient (r)

The Pearson’s correlation coefficient is a method that is based on using the covariance matrix of the data to evaluate the strength of the linear relationship between two variables. Normally, the Pearson’s correlation coefficient between two variables is [53]:

$$r\left(x,y\right)=\frac{{\sum }_{i=1}^n\left[\left(x_i-\overline{x}\right)\times \left(y_i-\overline{y}\right)\right]}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}\times \sqrt{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(1)

where $\overline{x}$ and $\overline{y}$ denote the means of the two variables, and n is the sample number (length of the time series).

In this study, the Pearson’s correlation coefficient was utilized to determine the response periods between SIF and climatic factors (P, VPD, ES, and RSM). For example, on each pixel over the grasslands, we first calculated the r between deseasonalized time series of SIF and P (r (SIF_ds, P_ds)). The deseasonalized P time series was moved one or more positions (one or more months) forward, then they were used to recalculate the correlation with deseasonalized SIF time series (r (SIF_ds, P_ds1), r (SIF_ds, P_ds2), …). We compared these r to ascertain the response period of SIF to P, more specifically, if the r (SIF_ds, P_ds2) was the maximum, the response period was between two and three months (i.e., the lag phase of SIF to P was two months).

In this study, all maps and figures were created by Python 3.7, and the cartopy (https://pypi.org/project/Cartopy/, accessed on 4 June 2021) and matplotlib (https://matplotlib.org/, accessed on 4 June 2021) libraries. The hardware and storage employed were the Inter(R) Xeon(R) Gold 5218 CPU @ 2.30 GHz 2.29 GHz and 192 GB Random Access Memory with 2 TB Read-Only Memory. The operating system was 64-bit Windows 10.

3. Results

3.1. The Impact of Seasonality on the Relationship of SIF to P, VPD, ES and RSM Dynamics

Firstly, we calculated the spatial r maps of SIF with P, VPD, ES, and RSM, without any time series adjustment (Figure 2). According to the mechanism relationship between SIF and the four climatic factors, P, ES, and RSM have positive correlations with SIF, and VPD has a negative correlation with SIF. However, in Figure 2, the spatial r of original SIF with the original P (r (SIF, P)), VPD (r (SIF, VPD)), ES (r (SIF, ES)), and RSM (r (SIF, RSM)) did not follow the mechanism. More specifically, in the spatial maps of r (SIF, P) (Figure 2a), r (SIF, ES) (Figure 2c), and r (SIF, RSM) (Figure 2d), there were large areas of negative r values, while in the spatial maps of r (SIF, VPD) (Figure 2b), the r values were almost positive.

We conducted further research and found that the seasonality/periodicity of the time series led to the anomalous phenomenon. We compared the original and deseasonlized time series between SIF and RSM at 50°N and 70°E (Figure 3) and found the original time series of SIF and RSM (Figure 3a) had clear seasonality and the dynamic response between SIF and RSM were inconspicuous, while Figure 3b showed a legible response between the deseasonlized time series of SIF (SIF_ds) and RSM (RSM_ds) dynamics. These indicated that seasonality affected the qualification of the relationship between SIF and the four factors, and even led to incorrect results. Therefore, for qualifying the relationship between SIF and the four climatic factors, the impacts of seasonality need to be removed. Simultaneously, the deseasonalization based on the STL method was a useful means of removing the disturbance generated by periodicity.

After deseasonalizing the time series of SIF, P, VPD, ES, and RSM (SIF_ds, P_ds, VPD_ds, ES_ds, and RSM_ds) based on the STL method, we recalculated the spatial correlation between SIF_ds and the four factors (Figure 4). Overall, the spatial r maps followed the mechanism mentioned above, and ES correlated with SIF best. However, there were still unreasonable r values between SIF and the four factors. For example, between the 90°E and 100°E, in the northern and southern regions of the grasslands, the r (SIF_ds, VPD_ds) were greater than zero, which they should had been negative correlations between SIF and VPD according to the response mechanism between them. In addition, the corresponding spatial maps of p showed a large area which did not pass the significance test (p > 0.05). The area with an insignificant r between P and SIF (Figure 4e) was the greatest, while that between ES and SIF (Figure 4g) was the lowest. The p value maps of SIF with VPD (Figure 4f) and RSM (Figure 4h) had analogous spatial distribution.

3.2. The Influence of Response Period on r of SIF with P, VPD, ES and RSM

Certainly, the delayed response of SIF to the four climatic factors means there are time lag effects between them. By using the methods mentioned in Section 2.3, we determined the lag phases of SIF to the four factors (Figure 5, masked by the significant extremum r of SIF to the four climatic factors). SIF had a time lag effect on the four factors over the grasslands and the lag phases were mainly within one month. Combing Figure 4 and Figure 5, the spatial r of SIF with P, VPD, ES, and RSM were consistent with the spatial lag phase maps of SIF to the four factors. Specifically, poor (p > 0.05) or opposite correlations between SIF and these factors were mainly located in areas with lag phases greater than one month. Collectively, the r values were affected by the response periods and the longer the period, the greater the impact. In addition, in Figure 5, the area with response periods of SIF to ES of more than one month were the lowest among the four factors, which meant the response of SIF to ES was timely.

We integrated the maximum r of SIF with P, ES, and RSM and the minimum r of SIF with VPD, and the corresponding spatial p maps were also integrated (Figure 6). Compared with Figure 4, the spatial r maps were more reasonable, and the qualification of the relationship between SIF and climatic factor was improved (by comparing the colors). The overall spatial r of SIF with P, VPD, ES, and RSM in Figure 6 showed improvement, especially in the areas with poor or even inverse correlations. Additionally, the areas with p > 0.05 were noticeably decreased. Since the extremum r values were obtained by adjusting the time series, this indicated that the delayed response of SIF to the four climatic factors, or the time lag between them, impacted the assessment and qualification of the relationship among variables.

In addition, we conducted comparative statistical analyses (shown in Figure 4 and Figure 6), and the results (Figure 7) showed that all of the r of SIF with the four factors were improved. The improvement of the r of SIF with P and ES were the largest and the smallest, respectively, which suggested the response period affected the correlation of SIF to P the most and the correlation to ES was the least affected among the four factors. The response period, or time lag effect, also has different degrees of influence on evaluating and qualifying the relationship between vegetation growth or ecosystem productivity and the climatic factors. All of these indicated that the correlation of SIF with P, VPD, ES, and RSM were impacted by the time lag effect. In addition, though we corrected the spatial r of SIF with the four variables, there were poor r values in some areas; the reason for this may be that vegetation greenness and productivity are affected by multiple variables, and a single variable cannot fully explain the vegetation greenness and productivity dynamics.

4. Discussion

4.1. Methods of Deseasonlization

Section 3.1 showed that the seasonality/periodicity seriously impacted the qualification of the relationship between SIF and the climatic factors. We eliminated the seasonality disturbance through the STL method. In fact, differencing and classical decomposition methods were also widely used to deseasonalize the seasonality of the time series [54], and these methods can be implemented by computer languages such as Python and R. However, both of them have weaknesses. Differencing refers to using the next cycle minus the time series of this cycle, which will reduce a cycle of the observations. For the data of this study, a cycle is a year, and the difference method will thus reduce the length of the original time series by one year. For classical decomposition methods, there are two forms of classical decomposition: an additive decomposition and a multiplicative decomposition [55]. They assume that the seasonal component repeats from year to year, which means they are unable to capture these seasonal changes over time. Classical decomposition methods are also unable to handle outliers properly. The STL method does not have these problems, and it can be robust to outliers. Although the STL method does not handle trading day or calendar variation automatically and only provides facilitation for additive decompositions, it does not affect the extraction of the seasonality of SIF and the four factors [52]. We believe that the STL method can play a greater role in the field of natural sciences.

4.2. The Determination of Lag Phase

In this study, by comparing the maximum or minimum value of Pearson’s correlation coefficient (r), we determined the lag phase, which is inspired by the studies of Peng et al. [56] and Xie et al. [57]. Our results showed a delayed response or time lag effect which impacted the correlation between SIF and the climatic factors. However, there were some uncertainties in this method. More specifically, Figure 6 could be considered as correcting the r of SIF with the four factors through adjusting the time series of SIF and the four factors to a reasonable response period. Compared with Figure 4, the spatial maps of p were noticeably improved in Figure 6. In fact, we adjusted the time series within 6 months, but in some areas the r of SIF with the climatic factors were insignificant at any adjustment. For this situation, we cannot determine the lag phases between variables, or even analyze the response between them. In addition, the correlation coefficient is a measure of the overall relationship between two variables. However, in the study period, some events, such as droughts and floods, made the relationship between SIF and the climatic factors abnormal or even deviant. This also affected the assessment and qualification of the mechanism of the relationship between variables. For determining a more accurate lag phase or response period, we believe that it is more effective to study the response between the variables from the perspective of events and utilize higher temporal and spatial resolution data.

4.3. The Impact of Time Lag Effect on Ecosystem and Climatic Change

The time lag effect of major climatic driving forces on phenology should be paid more attention to improve the prediction and evaluation of vegetation response [56,58,59]. An increasing number of studies have found that the responses of vegetation to climate has a certain time lag [60,61,62,63,64,65]. Braswell et al., researched the interannual lag of vegetation’s response to temperature among various ecosystems [30]. Wang et al., found the spatial pattern of the normalized difference vegetation index (NDVI) is strongly correlated with the longer time intervals of precipitation than with recent precipitation [66]. Anderson et al., investigated the time-lag effects for different climatic factors (i.e., radiation, rainfall, and aerosol optical depth) in the Amazon and found that the time lag of the vegetation responses to the climatic factors differed in the rainforest [58]. The vegetation responses to soil moisture lagged by approximately one month over the Australian mainland [64]. These studies recognize the importance of the time lag effects of vegetation responses to climate. The four climatic factors, P, VPD, ES, and RSM, were related to the vegetation water stress. Our findings showed the SIF mainly lagged behind these factors by one month, and in small areas the lag phases were within two and three months. This was similar to the results of Zhao et al. [67], which found the NDVI lagged behind standardized precipitation and the evapotranspiration index (SPEI, a drought index which can characterize water stress) within 2–3 month over grassland vegetation in the Chinese Loess Plateau (102°54′–114°33′E, 33°43′–41°16′N). The results of Zuo et al., displayed that the NDVI has a one-month lag effect on precipitation over the grassland of the Yarlung Zangbo River Basin (82°00′–97°07′E, 28°00′–31°16′N), which is the same as our results in this region [68]. Collectively, vegetation growth may not be primarily driven by present climatic conditions, but earlier climatic conditions may also have impact on vegetation growth. Therefore, the time lag effects should be considered when exploring the mechanisms of climate-vegetation interactive effects.

Actually, most studies that focused on large-scale climate-vegetation interactions adopted the simultaneous climatic and vegetation indicators, and the time lag effects were minimally considered, which would increase the uncertainty of the results. For example, the composite drought index involved in hydrological, meteorological, and ecological variables only combined these variables by methods such as linear combination, regression, and weight combination, and did not consider the time lag effects between them. When using drought indices to monitor and assess vegetation drought, most studies did not consider the response mechanism or time lag effect, leading to uncertainties in the results [26]. Additionally, for the estimation of terrestrial variables in ecological and hydrological process [69], for instance evapotranspiration, the estimation may involve more variables than were used in building a composite drought index. The time lag effects among various variables also need to be carefully considered in related models, which may improve the accuracy of the estimate. In fact, as climate change intensifies, climatic system become anomalous. We need more accurate results, whether in the monitoring and prediction of natural disasters or the estimation of surface variables. Taking the response mechanisms or time lag effect into account may be a potential direction to achieve this end.

4.4. Limitations and Prospects

There were limitations in this study. The length of the response period may not be accurate enough, since the temporal resolution of the datasets we used were monthly. With regard to the spatial r of SIF with the four variables, there were insignificant or poor r values in some areas after adjusting the time series of variables. The potential reason for this was that a single climatic variable change cannot fully explain the vegetation greenness or productivity dynamics. We did not comprehensively utilize P, VPD, ES, and RSM to explore and analyze the response of SIF to a combined indictor or index. In future studies, we will combine multiple climatic factors to generate a composite indictor to explore its relationship with vegetation growth, ecosystem productivity, or carbon exchange. We will also utilize higher temporal and spatial resolution datasets to determine the response period of the ecosystem to climatic factors. We will also extend the study area to larger regions with diverse vegetation types in order to analyze the dynamic response between different vegetation types and climatic factors.

5. Conclusions

In this study, we utilized SIF, P, VPD, ES, and RSM to explore the response of vegetation greenness or productivity to the climatic factor dynamics over the Eurasian Steppe. Our findings showed the seasonality of the variables’ time series created disturbance in qualifying the relationship between SIF and the climatic factors, and the deseasonalization based on the STL method was an effective means to remove the disturbance of periodicity. The response of SIF to P, VPD, ES, and RSM all had time lag effects, and the response period had an impact on the correlation of SIF with the four climatic factors. Among P, VPD, ES, and RSM, P was the most affected by the response period, while ES was the least impacted by the response period. This study provides a reference for researchers to attribute and monitor the dynamics of ecosystem productivity and carbon exchange.