Non-uniform seasonal warming regulates vegetation greening and atmospheric CO₂ amplification over northern lands

There are 19 CO₂ measurement sites in the NOAA's Global Greenhouse Gas Reference Network (https://esrl.noaa.gov/gmd/ccgg/ggrn.php) and 3 sites in Scripps CO₂ program (http://scrippsco2.ucsd.edu/data/atmospheric_co2/sampling_stations) located at lands over 50° N. Among these sites, only the Barrow (BRW) observatory site recorded the atmospheric CO₂ concentration (i.e. [CO₂]) continuously during 1982–2010. Thus, in this study, the in situ long-term CO₂ measurements from BRW were regarded as the homogeneous CO₂ concentration in northern lands (>50° N). The seasonal curve of this [CO₂] record was shown in figure S1 (available online at stacks.iop.org/ERL/13/124008/mmedia).

The in situ [CO₂] observations at the BRW site were collected hourly. The daily and monthly data provided by NOAA were averaged from the hourly observations. This study used the monthly data to derive the CO₂ amplitude. The anomalies of monthly [CO₂] (i.e. monthly [CO₂] − yearly mean [CO₂]) were first calculated and then were used to derive the maximum (i.e. [CO₂]_max) and minimum (i.e. [CO₂]_min) monthly [CO₂] in each year. The difference between [CO₂]_max and [CO₂]_min was defined as the CO₂ amplitude ([CO₂]_amplitude).

To avoid the biases of various processing algorithms, we collected the estimates of [CO₂]_amplitude from Forkel et al (2016) and the GLOBALVIEW products (figure S2). In Forkel et al (2016), the time series of daily [CO₂] records were first fitted with polynomial and harmonics functions and then de-trended with the Fast Fourier Transformation method (Thoning et al 1989, Thoning et al 2015). In each calendar year, the [CO₂]_amplitude was calculated as the peak-to-trough difference of the de-trended seasonal cycle. The data in the GLOBALVIEW-CO₂ product has already been smoothed, interpolated and extrapolated (Masane and Tans 1955). GLOBALVIEW-CO₂ provides observations at 7 d intervals. We obtained its [CO₂]_amplitude as the difference between the maximum and minimum weekly CO₂ data in each calendar year. According to the Theil–Sen analysis, the long-term trends of [CO₂]_amplitude were consistent among these processing methods (figure S2).

The temperature trends were analyzed based on the latest version of the CRU temperature data (CRU TS4.0). It is gridded with a spatial resolution of 0.5° × 0.5° at a monthly time step. This product is gridded using the Angular-distance weighting interpolation (Harris and Jones 2017) based on the observations collected from 2600 stations worldwide (Harris et al 2014). The CRU climate products have been widely used for phenology analysis and for driving different types of ecosystem models (Koven et al 2011, Fu et al 2014, Xia et al 2017). In this study, seasonal temperature was averaged from monthly temperature following the definition of four seasons: Spring, March–May; Summer, June–August; Autumn, September–November; Winter, December–February. Besides, MAT of a given year was defined as the average of the monthly temperature from January to December. The Theil–Sen estimator and Mann-Kendall trend test were applied in detecting the temporal trends of the seasonal temperature (see more details in section 2.6).

The normalized difference vegetation index (NDVI) is widely used as an indicator for vegetation productivity (Myneni et al 1997, Zhou et al 2001). It is calculated as the normalized ratio between near infrared and red reflectance bands (Tucker 1979, Tucker et al 2005). The NDVI used in this study is from the Advanced Very High Resolution Radiometer sensors, which has the longest record of continuous satellite data since 1981. Here, we used the newest version of GIMMS NDVI dataset (NDVI3g) (Tucker et al 2005, Pinzon and Tucker 2014). It is a global product at spatial resolution of ∼8 × 8 km and temporal resolution of 15 d. The NDVI3g has been widely used for analyzing vegetation changes in recent years (Tucker et al 2005, Peng et al 2013, Wang et al 2014). The maximum value composites method (Holben 1986) was applied to merge segmented data strips to half-monthly values. To lessen the impacts of sparse soils and snows on vegetation, following Zhang et al (2013), the areas with multiyear average NDVI less than 0.1 in the northern lands (>50° N) were removed from the analysis. Moreover, only the half-monthly NDVI from January to September were used to derive the phenology indices (i.e. start and end of growing season length (GSL), see section 2.4).

The sum of monthly NDVI from January to December in a certain year was regard as the yearly NDVI (Gonsamo et al 2016, 2017). Note that the illegitimate signals of the half-monthly NDVI data (i.e. NDVI < 0.1) were filtered. Regional NDVI used in the trend analysis were averaged from the grids of the northern lands (>50° N). Before calculating the sensitivity of NDVI to MAT and the partial correlation between NDVI and seasonal temperature (see section 2.6), yearly NDVI data were resampled to raster of 0.5° × 0.5°, to couple with the CRU temperature data.

GSL was calculated as the difference between the start (SOS) and end (EOS) of growing season. The SOS and EOS were retrieved from the seasonal NDVI curve in each year based on the NDVI green-up thresholds, which were determined from the rate of seasonal changes in the multiyear mean NDVI (Piao et al 2006, 2011, Zhang et al 2013). More specifically, there were six steps in the determination of GSL. First, we calculated the seasonal curve of multiyear mean NDVI from 1982 to 2010 for each land grid cell and obtained the changing rate of NDVI (NDVI_ratio) as:

$$\begin{eqnarray*}&&\begin{array}{l}{{\rm{NDVI}}}_{{\rm{ratio}}}(t)=[{\rm{NDVI}}(t+1)-{\rm{NDVI}}(t)]\\ \,\,\,\,\,\,/[{\rm{NDVI}}(t)],\end{array}\end{eqnarray*}$$

where t is time throughout the year with an interval of 15 d. Then, after removing evident noise in the multiyear mean time-series curve of NDVI for each land grid cell, we performed a least-square regression analysis on the curves from January to September and from July to December for determining the NDVI thresholds of SOS and EOS, respectively, with an inverted parabola equation:

$$\begin{eqnarray*}&&{\rm{NDVI}}=a+{a}_{1}x+{a}_{2}{x}^{2}\,+\,\ldots \,+\,{a}_{n}{x}^{n},\end{eqnarray*}$$

where x is the Julian days and n is 6. The corresponding NDVI(t) with the maximum or minimum NDVI_ratio was used as the NDVI threshold for determining SOS or EOS, respectively. Next, we performed a least-square regression analysis on the NDVI time-series curve in each year for each pixel after removing noise in the two different periods. After that, the SOS and EOS were identified from the fitted NDVI seasonal curves and their NDVI thresholds, by selecting the day when the fitted NDVI curve first reached the NDVI threshold. Finally, the GSL in each year for each land cell was calculated from the difference between EOS and SOS. The method in this study had been validated by ground-based phenological data in Tibetan Plateau (Zhang et al 2013) and has been widely used for detecting phenological changes in various regions (Piao et al 2006, Piao et al 2011, Zhang et al 2013).

Outputs of annual GPP from five TBMs which provided all the land cells above the 50° N in the model integration group of the Permafrost Carbon Network (http://permafrostcarbon.org/) were analyzed in this study. The five TBMs are UVic (Peter 2001, Matthews et al 2004), CoLM (Dai et al 2003, Ji et al 2014), CLM4.5 (Oleson et al 2003), TEM604 (Hayes et al 2011) and ORCHIDEE (Krinner et al 2005). Details about these TBMs were listed in table S1. The simulation protocol and model's driving data have been described in previous studies (Rawlins et al 2015, McGuire et al 2016, Peng et al 2016). A flux-tower-based GPP database was also used in this study. It was up-scaled from FLUXNET observations (44 sites locating in the lands northern 50° N) of carbon dioxide, energy and water fluxes with the machine learning technique of model tree ensembles (MTE, Jung et al 2011). The MTE GPP is a global gridded product with a resolution of 0.5° × 0.5°. This product has been widely used as benchmarks to evaluate model performance in recent years (Anav et al 2013, Tjiputra et al 2013, Peng et al 2015, Xia et al 2017).

We estimated the linear trends of the CO₂ indices (i.e. [CO₂]_amplitude, [CO₂]_max and [CO₂]_min), vegetation and temperature using a non-parametric Theil–Sen estimator over each time period. The significance of the trend was computed by the Mann-Kendall trend test. Comparing with the ordinary least squares estimation, the Theil–Sen estimator and Mann-Kendall trend test is less sensitive to outliers (Fernandes and Leblanc 2005, Wang et al 2018). The temporal anomalies were used for the linear-trend analyses. This analysis can provide both trend and its level of significance (i.e. the P value that quantifies the probability of whether the trend is statistically significant from zero) for each period.

The moving-window method was used to detect whether the increasing trends of CO₂ indices are persistent. Comparing with the piecewise linear fitting method, it less depends on the results of single linear segment and the interval-length (Schleip et al 2008). This method has been used in detecting the changes of growing-season length and its response to climate change on various time scales (Rutishauser et al 2007, Schleip et al 2008, Jeong et al 2011, Fu et al 2015). Because the results based the moving-window analysis may be affected by the window-length (Fu et al 2015), we repeated the moving-window analyses with ten-, 15- and 20-year lengths.

The temporal trends of MAT (${{\rm{\Delta }}}_{{\rm{M}}{\rm{A}}{\rm{T}}}$), the sensitivity of NDVI to temperature (${\gamma }_{{\rm{N}}{\rm{D}}{\rm{V}}{\rm{I}}}$), and the sensitivity of CO₂ amplitude to temperature (${\gamma }_{{[{{\rm{CO}}}_{2}]}_{{\rm{a}}{\rm{m}}{\rm{p}}{\rm{l}}{\rm{i}}{\rm{t}}{\rm{u}}{\rm{d}}{\rm{e}}}}$) in the periods of 1982–2010 and 1993–2007 were calculated in each grid cell. The ${\gamma }_{{\rm{N}}{\rm{D}}{\rm{V}}{\rm{I}}}$ and ${\gamma }_{{[{{\rm{CO}}}_{2}]}_{{\rm{a}}{\rm{m}}{\rm{p}}{\rm{l}}{\rm{i}}{\rm{t}}{\rm{u}}{\rm{d}}{\rm{e}}}}$ were derived from the slope of linear regressions, representing the changes of NDVI and [CO₂]_amplitude with per degree change of MAT. One-way ANOVA was used to estimate their differences between 1982–2010 and 1993–2007. All the analyses were applied in R (http://r-project.org/).

Partial correlation analyses were applied to exclude the impacts of the co-varying factors. For example, in calculating the impact of spring temperature on SOS, spring precipitation, spring solar radiation and last-year's autumn temperature were set as the controlling variables. Autumn precipitation, autumn solar radiation and spring temperature were set as the controlling variables to quantify the impact of autumn temperature on EOS. Similar method was used to detect the impact of spring and autumn temperature on annual NDVI and modeled GPP by replacing the seasonal precipitation and solar radiation with annual values. All the data were aggregated to the 0.5°× 0.5° resolution. The precipitation data were derived from the CRU TS4.0 dataset (Harris et al 2014) and the radiation data was from the Terrestrial Hydrology Research Group at Princeton University (Sheffield et al 2006).

We first examined the trends of the measured annual CO₂ amplitude ([CO₂]_amplitude) at Point BRW, Alaska (71°N). The non-parametric Theil–Sen estimator showed that the increasing trends of the [CO₂]_amplitude at the BRW (0.075 ppm yr⁻¹, P < 0.05; figure 1(a)) were associated with the decreasing [CO₂]_min (−0.058 ppm yr⁻¹, P < 0.05) rather than the enhancing [CO₂]_max (0.016 ppm yr⁻¹, P = 0.17) from 1982 to 2010. The ten-year moving windows show that the increasing rates of the CO₂ amplitude was slower around 2000 (figure 1(a) and table S2). To avoid the biases from different time-intervals for trend estimation, we also detected the trends with 15-year (figure S3) and 20-year (figure S4) moving windows. The results also showed that the trends of [CO₂]_amplitude from mid-1990 to mid-2000 (e.g. 0.03 ppm yr⁻¹, P = 0.55, 1993–2007) (figure S4(a) and table S4) were significantly slower than those during other periods. A recent study which integrated the CO₂ records from multiple sites also has showed a stalled trend in the seasonality of atmospheric CO₂ during the same period (Yuan et al 2018).

Zoom In Zoom Out Reset image size

A slowdown of vegetation greening since mid-1990s was also observed by our analysis on the dynamics of NDVI (figures 1(b) and S5(d)). This finding is consistent with the results from recent analyses on vegetation dynamics over the pan-Arctic tundra (Bhatt et al 2013, Jeong et al 2013). Both the MTE GPP (figure S5(e)) and the ground-based measurements of growing-season net ecosystem CO₂ exchange (figure S5(f)), the net ecosystem exchange data were derived from Belshe et al (2013) showed similar trends since 1990s. These lines of evidences together suggest that the increasing trend of the growing-season CO₂ uptake has weakened from 1990s to mid-2000s in northern ecosystems.

Meanwhile, a weak but significant linear relationship between average NDVI over the northern lands (>50° N) and [CO₂]_min at the BRW site was observed during 1982–2010 (r = −0.47, P < 0.05; figure 1(b)). The synchronous changes of [CO₂]_amplitude with [CO₂]_min (figures S5(a)–(c)) and NDVI (figure 1) imply that the long-term positive trend of [CO₂]_amplitude is, at least in part, driven by photosynthetic CO₂ uptake or vegetation growth (Forkel et al 2016, Wenzel et al 2016).

An increasing body of research has shown a non-linear response or reduced sensitivity of vegetation growth to rising MAT over high latitudes in recent years (Bhatt et al 2013, Jeong et al 2013, Piao et al 2014). As shown by figure 2, although the MAT increased even faster from mid-1990s to mid-2000s (e.g. 1993–2007) than the whole period of 1982–2010, the sensitivities of [CO₂]_amplitude and NDVI to MAT were lower during that period than 1982–2010. Given the fact that temperature in different seasons has non-uniform impacts on vegetation growth (Xia et al 2014), we further analyzed the changes of seasonal temperatures based on the CRU temperature datasets (see methods). As shown by figure 3(b), the fastest warming season was spring during 1985–1999 (+0.12 °C year⁻¹) but then changed to autumn during mid-1990s to mid-2000s (e.g. +0.11 °C year⁻¹ in 1993–2007, table S4). It indicates that a better understanding of the relationship between the seasonal temperature changes and vegetation growth is needed to explain the slowdown of [CO₂]_amplitude from 1990s to mid-2000s.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The variation of NVDI during 1982–2010 in northern ecosystems depends substantially on the GSL on both grid and regional scales (figures 1(b) and S6). Further partial correlation analysis showed that the SOS (partial r = −0.36) was more dependent on temperature change than the EOS (partial r = 0.018, figures 3(c) and S7). During 1982–2010, the SOS was advanced by 2.15 d °C⁻¹ with spring warming, whereas warming in autumn only delayed EOS by 0.80 d °C⁻¹ in northern ecosystems. As a result, the advancing rate of SOS over the moving 15 years has decreased but the extending rate of EOS was not significantly increased during 1982–2010 (figure S8 and table S5). Meanwhile, the increasieng rate of NDVI until it stalled in mid-1990s is driven by warming-induced increase in spring and early summer NDVI along with the advancement of SOS (figure S8), for which spring warming has stalled after mid-1990s. It has contributed to the decline in the rate of [CO₂]_amplitude increase since the mid-1990s. These results together suggest that the non-uniform warming between spring and autumn during mid-1990s and mid-2000s (figure 3(b)) could be an important driving factor for the slowdown of expanding GSL, greening vegetation and the decreasing [CO₂]_min. This can qualitatively explain the pause in the enhancement of [CO₂]_amplitude (figure 1(a)).

We examined whether TBMs that have been focusing on simulation of C dynamics in northern latitudes can adequately represent the differential impacts of spring and autumn warming on vegetation productivity. The ensemble output of GPP from five TBMs (CLM4.5, CoLM, ORCHIDEE, TEM6 and UVic; table S1) and a flux-based GPP dataset (MTE) were analyzed (figure 4). The dependence of modeled GPP variations on spring-temperature change (with the inter-model mean partial r as 0.57 under the significant level of P < 0.05) is comparable with that of the MTE GPP (partial r = 0.53, P < 0.05) as well as that of NDVI (partial r = 0.49, P < 0.05; figure 4(b) and figure S9). However, the dependences of GPP variations on autumn-temperature change is more positive in the models (inter-model mean partial r = 0.31, P < 0.05) than the MTE GPP (partial r = −0.04, P < 0.05) and NDVI (partial r = −0.16, P < 0.05; figures 4(b) and S10). This mismatch between modeling and data-oriented results indicates that the current TBMs overestimate the positive impact of rising MAT on ecosystem CO₂ uptake in the autumn.

Zoom In Zoom Out Reset image size

As the land biogeochemical component in most Earth system models is similar to the TBMs in this study, it is still challenging to accurately simulate the seasonal cycle of atmospheric CO₂. The findings of this study suggest that a better representation of the warming impacts on autumn phenology could partially improve the models' performance. However, the autumn phenology is diversely represented in different models. For example, leaf senescence in the ORCHIDEE model is simulated as the timing when monthly temperature falls below a given number, which varies with plant function type (Krinner et al 2005). In the TEM, growing season ended when the soil temperature is lower than the frozen point. However, leaf-senescence events are collectively affected by not only temperature but also day length (Ballantyne et al 2017), radiation (Bauerle et al 2012) and even spring phenology (Keenan and Richardson 2015, Liu et al 2016). In fact, the poor representation of autumn phenology by the models has been raised in some previous studies (Richardson et al 2010, 2012, Keenan and Richardson 2015). Thus, combining the different types of phenological data (e.g. Richardson et al 2018) with better phenology models could be helpful to improve the simulation of the seasonal cycle of atmospheric CO₂ at high latitudes in Earth system models.

This study highlights that the asymmetric responses of vegetation growth to spring and autumn warming is an important driver for the decadal changes in the seasonality of atmospheric CO₂. In spring, solar radiation is not limiting as temperature (Tanja et al 2003). Thus, spring warming extends the GSL by advancing the onset of plant photosynthesis (Piao et al 2008), leading to the increasing vegetation productivity and the decreasing [CO₂]_min of the atmospheric CO₂ seasonal cycle. In autumn, solar radiation can obstruct the accumulation of abscisic acid (Gepstein and Thimann 1980) and substantially delay the timing of leaf senescence. Photoperiod is a more proximal factor than temperature in controlling senescence (Bauerle et al 2012). Thus, autumn warming has a neutral impact on vegetation productivity (figure 3(c)) over the northern lands.

Warming in autumn as well as in spring could potentially enhance the peak of atmospheric CO₂ seasonal cycle by stimulating the respiratory processes of plants and soil microorganisms (Piao et al 2008, 2017a). As shown by the FLUXCOM database (Jung et al 2017), the increasing trends of total ecosystem respiration during both growing and non-growing seasons were significantly larger in mid-1990s to mid-2000s (e.g. 1993–2007) than 1982–2010 (figure S11). This result is consistent with previous findings that the warming induced increases in respiration could partially cancel out the impact of enhanced photosynthesis on the atmospheric CO₂ seasonality in North Hemisphere (Gonsamo et al 2017, Jeong et al 2018). However, further conclusions are limited by quantifying the contributions of increased respiration to the slowdown of CO₂ amplitude since mid-1990s. Future studies could improve on the present analysis through breaking the limitation.

Both spring and autumn are likely to keep warming in future scenarios (IPCC 2013), and further warming could trigger some limitations on vegetation productivity. For example, early spring warming may slow the fulfillment of chilling requirement for spring leaf phenology and thus delay the SOS (Yu et al 2010, Fu et al 2015, Vitasse et al 2018). The spring warming induced advancement of leaf unfolding date could increase the risk of frost damage to buds (Inouye 2008) and decrease soil water availability for subsequent peak productivity (Buermann et al 2013). Autumn warming may cause more cloudy days with less radiation (Vesala et al 2010), which may accelerate the ending of growing season (Bauerle et al 2012). Meanwhile, warm autumns strengthen the evapotranspiration during the late growing season and intensify the stresses of drought on vegetation growth (Barichivich et al 2013).