Increasing terrestrial ecosystem carbon release in response to autumn cooling and warming

Climate data

Monthly climate data (air temperature at 2 m and cloudiness) with a spatial resolution of 0.5° were obtained from the CRU Time Series 4.0.¹⁵ We extracted data from 1982 to 2018 to match the time series of satellite vegetation observations. The VPD was calculated as the difference between saturated water-vapour pressure and actual water-vapour pressure³¹. Temperature and vapour-pressure data used for the VPD calculation were obtained from CRU.

Soil moisture data

The daily root-zone soil moisture with a spatial resolution of 0.25° for the period 1980–2018 was obtained from the Global Land Evaporation Amsterdam Model (GLEAM v.3.3a)³². The dataset is based on radiation and air temperature from a reanalysis, a combination of gauge-based, reanalysis-based and satellite-based precipitation and satellite-based vegetation optical depth.

Fire emission data

Monthly carbon emissions from biomass burning were obtained from the fourth-generation Global Fire Emission Database³³. This dataset has a spatial resolution of 0.25° and provides global data on the burning area and emissions on three-hourly, daily and monthly timescales and estimates the contributions of different fire types. Emissions data can be obtained for different substances, such as carbon (C), dry matter (DM), carbon dioxide (CO₂), carbon monoxide (CO) and methane (CH₄).

Satellite vegetation greenness data

The satellite-based NDVI archived from the MODIS NDVI dataset with a spatial resolution of 0.5° and a temporal resolution of 16 days was used here to detect vegetation greenness changes. In addition, the solar-induced chlorophyll fluorescence product was used as a proxy of vegetation photosynthesis. We furthermore used the four-day clear-sky CSIF time series (2000–2019) with a spatial resolution of 0.05° × 0.05° from ref. ³⁴ (https://osf.io/8xqy6/).

GPP based on NIRv

The NIRv is a newly developed satellite vegetation index combining NDVI and near-infrared band reflectivity of vegetation and is recognized as a proxy of GPP^35,36. We obtained the 0.05° NIRv_GPP from 1982 to 2018 from ref. ³⁷. This product was produced by upscaling the relationships between NIRv and observed GPP to the global scale and was judged to perform well in capturing interannual trends of GPP³⁷.

Atmospheric CO₂ data

In situ observations of daily CO₂ concentration at Point Barrow were obtained from the National Oceanic and Atmospheric Administration/Earth System Research Laboratory network. According to analyses of atmospheric transport and mixing processes, the CO₂ signals detected at Barrow are suggested to be an integrated measure of carbon fluxes over both the high latitudes and the middle latitudes²⁰.

Ecosystem carbon fluxes

Simulations of ecosystem carbon fluxes (GPP, TER and NEE) derived from process-based model simulations (TRENDY), empirical models based on flux tower observations (FLUXCOM) and atmospheric CO₂ inversion models were jointly used for the investigation of net ecosystem carbon exchange over the northern middle and high latitudes.

The TRENDY dataset is an ensemble of dynamic global vegetation model (DGVM) simulations that are forced by CRU–National Centers for Environmental Prediction historical climate and CO₂ inputs³⁸. The DGVMs use a bottom‐up approach to simulate terrestrial CO₂ fluxes (for example, GPP, TER and NEE), and were extensively used to explore the mechanisms driving changes in carbon uptake and fluxes. The simulated GPP, TER and NEE from nine models of TRENDYv.8 (Supplementary Table 1) were used in this study. The S2 experiment, which considered the effect of both observed changes of CO₂ and climate on ecosystem carbon fluxes, was selected for studying the changes of ecosystem carbon fluxes before and after the temperature shift.

The FLUXCOM dataset is an upscaling product using empirical models forced by eddy-covariance data from 224 flux towers, remote sensing data and climate data^8,9,10. It provides estimates of global energy and carbon fluxes (http://www.fluxcom.org/). The empirical models were trained by three machine learning algorithms, including Random Forests, Artificial Neural Networks and Multivariate Adaptive Regression Spline, and thus provide a series of estimates of global carbon fluxes. We used the FLUXCOM carbon fluxes data driven by the European Centre for Medium-Range Weather Forecasts Reanalysis v.5 (ERA5) climate reanalysis from 1979 to 2018.

The atmospheric CO₂ inversion datasets provide estimates of NEE over land from long-term atmospheric CO₂ measurements using atmospheric transport models. Three atmospheric CO₂ inversion products were used here: monthly net biome production with a spatial resolution of 3.75° × 2.5° from the JENA CarboScope (version s76_vo2020) for the period 1976–2019, long-term global CO₂ fluxes estimated by the NICAM-based Inverse Simulation for Monitoring CO₂ (NISMON-CO₂) between 1990 and 2019 and the Copernicus Atmosphere Monitoring Service¹² (CAMS v.19r1) dataset between 1979 and 2019.

Eddy-covariance CO₂ observation data

The eddy-covariance measurements of carbon fluxes from tower sites were obtained from the Integrated Carbon Observation System 2018 and the FLUXNET Network 2015. We selected 48 eddy-covariance CO₂ observation sites with 10 yr continuous data (Supplementary Table 2) located north of 25° N and extracted temperature and NEE data from September to November to explore the change of ecosystem carbon exchange in autumn.

NEE estimation

The monthly NEE was estimated as the difference between TER and GPP. The autumn (September to November) GPP and TER derived from TRENDY and FLUXCOM over the study region were obtained by aggregating GPP and TER from each grid cell weighted by the grid-cell area. The NEE derived from atmospheric CO₂ inversions was directly used and compared against those from TRENDY and FLUXCOM. To compare the NEE before and after the temperature turning point, we divided the NEE time series into two periods: 1982–2003 and 2004–2018.

Calculation of the AZC

We used observations of CO₂ from Point Barrow to characterize the trends in the zero-crossing date of CO₂ (downward in spring and upward in autumn). These trends roughly correspond to the beginning of net carbon uptake in spring and the beginning of net carbon release in autumn. According to the method of ref. ³⁹, we obtained the detrended seasonal CO₂ curve by separating the seasonal cycle from the long-term trend and short-term variations, fitting a function consisting of a quadratic polynomial for the long-term trend and four harmonics for the annual cycle to the daily data. The residuals from this function fit are then obtained. A 1.5-month and a 390-day full-width half-maximum-value averaging filter were used for the digital filtering of residuals to remove the short-term variations and the long-term trend, respectively. Then we got the zero-crossing dates when the detrended seasonal CO₂ curve crosses the zero line from positive to negative and negative to positive, respectively.

The autumn carbon release is calculated as the amount of CO₂ released between the autumn zero-crossing date and the first week of September following ref. ²¹.

Identification of turning point of temperature

We used the piecewise linear regression method to determine the turning point of the mean autumn (September to November) temperature during 1982–2018 over the area north of 25° N. In addition, a moving t-test method was used to verify the turning-point identification. Then, the temporal trends of the mean autumn temperature before and after the turning point were calculated using the Mann–Kendall non-parametric trend test method, and the confidence intervals were determined using Sen’s slope statistics. According to the temperature trends before and after the turning point, we further identified the CAs as where the autumn temperature shows a decreasing trend after the turning point (2004) relative to that before the turning point, and WAs as regions outside the CAs. To maintain spatial integrity and continuity, we ignored the significance of the temperature trend when dividing the CAs and WAs.

To verify that our analysis is not affected by the division of the time period and regions, we also identified the temperature turning point at each grid point using the piecewise linear regression method and then extracted those grid points with significant temperature change and significant NEE change (P < 0.1) before and after the turning point. In addition, to ensure an adequate length of the time series for the trend analysis, only grid points with at least 15 yr of time series of temperature records and NEE data before and after the turning point were selected. For NEE, only TRENDY and CAMS had long enough time series to satisfy this criterion. Then, the selected grid points were classified into two groups, the warming group (shift from significant warming to a larger significant warming rate after the turning point) and the cooling group (shift from significant warming to significant cooling after the turning point), and the NEE trends before and after temperature turning point were compared.

Statistical analysis

Trend statistics and slopes of other climatic variables, carbon fluxes (GPP, TER, and NEE) and satellite indices during autumn were also determined using Mann–Kendall non-parametric trend test and Sen’s slope statistics. The Student’s t test was used to determine the significance of the trend difference in carbon fluxes during different periods.

The EOF decomposition method was used to explore the temporal and spatial variation of autumn temperature change. This method can determine the main components of variables on a time–space scale and provide the variance contribution of different distribution characteristics⁴⁰. Here, we decomposed autumn temperatures from 1982 to 2018 using EOF and selected the top three distribution modes according to the variance contribution degree, which represents the main spatial pattern and temporal changes of temperature from 1982 to 2018.

The actual value of NEE differs greatly among different sources of simulations. Consistent with previous studies, different sources of NEEs were normalized using the Z-score method to highlight the interannual variability and trend of NEE change.

Pearson’s correlation analysis and partial correlation analysis were performed to examine the relationship between climatic variables and carbon fluxes (GPP, TER and NEE) and satellite indices. We detrended all the variables before the partial correlation analysis using a linear regression method. The significance of correlations was assessed at P < 0.1.

We used a 15 yr sliding window to examine the dependence of the GPP–temperature relationship on mean autumn temperature. In each 15 yr window, the mean autumn temperature and the partial correlation coefficient between GPP and temperature were calculated. A total of 23 GPP–temperature correlation coefficients and the same number of mean autumn temperature values were obtained during 1982–2018, and then a linear regression equation was used to fit them. The dependence of the GPP–SM relationship on mean autumn temperature was examined using the same method.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com