Contrasting effects of CO₂ fertilization, land-use change and warming on seasonal amplitude of Northern Hemisphere CO₂ exchange

Continuous atmospheric CO₂ monitoring data indicate an increase in the amplitude of seasonal CO₂-cycle exchange (SCA_NBP) in northern high latitudes. The major drivers of enhanced SCA_NBP remain unclear and intensely debated, with land-use change, CO₂ fertilization and warming being identified as likely contributors. We integrated CO₂-flux data from two atmospheric inversions (consistent with atmospheric records) and from 11 state-of-the-art land-surface models (LSMs) to evaluate the relative importance of individual contributors to trends and drivers of the SCA_NBP of CO₂ fluxes for 1980–2015. The LSMs generally reproduce the latitudinal increase in SCA_NBP trends within the inversions range. Inversions and LSMs attribute SCA_NBP increase to boreal Asia and Europe due to enhanced vegetation productivity (in LSMs) and point to contrasting effects of CO₂ fertilization (positive) and warming (negative) on SCA_NBP. Our results do not support land-use change as a key contributor to the increase in SCA_NBP. The sensitivity of simulated microbial respiration to temperature in LSMs explained biases in SCA_NBP trends, which suggests that SCA_NBP could help to constrain model turnover times.

1 Introduction

The increase in the amplitude of seasonal atmospheric CO₂ concentrations at northern high latitudes is one of the most intriguing patterns of change in the global carbon (C) cycle. The seasonal-cycle amplitude (SCA) of atmospheric CO₂ in the lower troposphere at the high-latitude monitoring site of Point Barrow, Alaska, has increased by about 50 % since the 1960s (Keeling et al., 1996; Dargaville et al., 2002). Increasing SCA has also been registered at other high-latitude sites, mostly above 50^∘ N (Piao et al., 2017), and appears to be driven primarily by changes in seasonal growth dynamics of terrestrial ecosystems (i.e. net biome productivity – NBP), but uncertainty remains about the relative contributions from different continents and mechanisms.

Some studies proposed that the trend in SCA is primarily driven by increased natural vegetation growth and forest expansion at high latitudes due to CO₂ fertilization and climate change (Graven et al., 2013; Forkel et al., 2016; Piao et al., 2017). Others (Gray et al., 2014; Zeng et al., 2014) suggested that agricultural expansion and intensification resulted in increased productivity and thus enhanced the seasonal exchange in cultivated areas at mid-latitudes. However, evidence suggests that crop productivity stagnated after the 1980s in many regions in the Northern Hemisphere (Grassini et al., 2013), which is not reflected in SCA trends in recent decades (Yin et al., 2018).

Studies using land-surface models (LSMs) to attribute trends to the suggested processes usually convert simulated fluxes to CO₂ concentrations using atmospheric transport models (ATMs) and compare the results to in situ measurements (Dargaville et al., 2002; Forkel et al., 2016; Piao et al., 2017) or over latitudinal transects (Graven et al., 2013; Thomas et al., 2016). These studies have shown that LSMs systematically underestimated SCA trends, but it is not clear whether these biases are due to LSM uncertainties or due to trends or errors in the ATM (Dargaville et al., 2002). Piao et al. (2017) addressed these problems by designing systematic model experiments to compare observed CO₂ concentrations at multiple sites with ATM simulations forced by an ensemble of NBP from different LSMs and an ocean biogeochemistry model. Point Barrow was the only site where nearly all models accurately described the trend in SCA, while at other sites, LSMs generally captured the sign of the trend in SCA but either underestimated or overestimated its magnitude. Piao et al. (2017) further reported that CO₂ fertilization and climate change drove the increase in SCA for sites >50^∘ N, but at mid-latitude sites, land use, oceanic fluxes, fossil-fuel emissions and trends in atmospheric transport may have contributed to the SCA trends.

Attributing changes in the seasonal amplitude of atmospheric CO₂ to specific processes requires analysing net surface fluxes as a function of changes in gross fluxes (photosynthesis, respiration and disturbance), which LSMs can provide. However, quantifying a bias in CO₂ concentration at a given site from a bias in fluxes simulated by a land-surface model (LSM) is difficult, since the biases can be affected by many other factors such as transport model characteristics, forcing data used, etc. Atmospheric inversions provide a consistent framework for assimilating in situ CO₂ concentration observations to estimate net CO₂ surface fluxes while accounting for errors in the prior fluxes and for some errors in the ATM (Peylin et al., 2013). At large spatial scales, the trends in SCA can be related to trends in the seasonal amplitude of CO₂ fluxes (i.e. SCA of NBP – SCA_NBP). Such an approach has been used to analyse trends in net CO₂ uptake in boreal regions (Welp et al., 2016).

The spatiotemporal distribution of terrestrial and oceanic surface fluxes estimated by inversions provides thus direct insight about the regional patterns of SCA_NBP that is fully consistent with the amplitude of CO₂ concentrations in all stations of the observational network used and constitutes a direct benchmark for SCA_NBP simulated by LSMs. Here, we use top-down (inversions; Chevallier et al., 2010; Rödenbeck, 2005) and bottom-up (TRENDYv6 LSMs; Le Quéré et al., 2018) estimates of terrestrial CO₂ fluxes at northern extra-tropical latitudes between 1980 and 2015 to (i) assess the ability of those LSMs to simulate inversion-based trends in SCA_NBP, (ii) attribute the trends in SCA_NBP to specific regions in the Northern Hemisphere and (iii) attribute the relative importance of drivers using the ensemble model framework.

Trends in SCA_NBP from the inversions are based on multiple in situ measurements and therefore provide a reference (and respective uncertainty) for evaluating the regional attribution by LSMs. Regarding process attribution, LSMs allow separating the contribution of different drivers through factorial simulations. However, the attribution by LSMs cannot be easily validated, which is especially problematic given that LSMs underestimate the trends in SCA at the latitudinal scale (Thomas et al., 2016). Thus, we compare (i) the process attribution by LSMs (as e.g. in Thomas et al., 2016; Piao et al., 2017), (ii) the statistical attribution based on inversion fluxes, (iii) the statistical attribution based on LSMs fluxes, directly comparable to the inversion results, and (iv) the statistical attribution based on the differences between factorial simulations (cross evaluation of i and iii).

Our approach allows us thus to constrain SCA_NBP trends at hemispheric and regional scales from both top-down (inversions) and bottom-up (LSMs) methods and to evaluate the process attribution by LSMs using top-down estimates of SCA_NBP.

2 Data

2.1 Atmospheric inversions

The inversion of a transport model to infer surface fluxes from concentration measurements is an ill-posed problem due to the dispersive nature of transport in the atmosphere and to the finite number of available measurements. This ill-posedness can be compensated for by using some prior information about the fluxes to be inferred. This prior information also drives the separation between natural and fossil fuel emissions in the estimation. In order to illustrate the diversity of the inversion results, we take the example of two inversion systems that provide results for the study period between 1980 and 2015. We analysed monthly surface CO₂ fluxes estimated by the inversion systems from the Copernicus Atmosphere Monitoring Service (CAMS; Chevallier et al., 2005, 2010) and from Jena CarboScope (Rödenbeck et al., 2003; Rödenbeck, 2005). The two inversions used here solve for fluxes on their ATM grid, thus minimizing aggregation errors for large regions (Kaminski and Heimann, 2001).

The CAMS version r16v1 (http://atmosphere.copernicus.eu/, last access: 16 November 2017) (Chevallier, 2017) provides estimates of ocean and terrestrial fluxes at 1.9^∘ latitude by 3.75^∘ longitude resolution. The CAMS inversion system assimilates observations from a variable number of atmospheric CO₂ monitoring sites (119 in total, providing at least 5 years of measurements) and uses the transport model from the LMDz general circulation model (LMDz5A) nudged to ECMWF-analysed winds. More details can be found in Chevallier et al. (2010).

The CarboScope v4.1 (available at http://www.bgc-jena.mpg.de/CarboScope/?ID=s, last access: 26 March 2018) provides several versions that assimilate a temporally consistent set of observations. We used these versions for the study period (1980–2015) to test the influence of the number of assimilated sites on the results. The s76, s85 and s93 versions have assimilated observations from 10, 23 and 38 sites since 1976, 1985 and 1993, respectively. Surface fluxes (ocean and land) are provided at the latitude–longitude resolution of 4^∘ × 5^∘ of the TM3 atmospheric transport model is used (Rödenbeck, 2005). In this version, the atmospheric model is forced by the National Centers for Environmental Prediction (NCEP) meteorological fields.

CarboScope further provides a sensitivity analysis of the s85 version fluxes to different parameters of the inversion. The sensitivity tests performed are as follows: “oc” – fixing the ocean prior, “eraI” – forcing the inversion with fields from ERA-Interim reanalysis instead of NCEP, “loose” and “tight” – scaling the a priori sigma for the non-seasonal land and ocean flux components by 4 (dampening) and 0.25 (amplification), respectively, “fast” – reducing the length of a priori temporal correlations, and “short” – reducing the length of a priori spatial correlations. The resulting latitudinally integrated SCA_NBP and respective trends are shown in Fig. S2 in the Supplement.

Since CAMS includes a larger, but time-varying, number of multi-year air-sampling sites as they are available, it constrains better spatial patterns, while CarboScope keeps a fixed set of sites covering a given period, using fewer sites but avoiding artefacts in the time series related to the appearance or disappearance of measurement sites.

2.2 Land-surface models

Land-surface models (LSMs) provide a bottom-up approach for evaluating terrestrial CO₂ fluxes (i.e. NBP) and allow deeper insight into the mechanisms driving changes in C stocks and fluxes. The TRENDY intercomparison project compiles simulations from state-of-the-art LSMs to evaluate terrestrial energy, water and CO₂ exchanges starting from the pre-industrial period (Sitch et al., 2015; Le Quéré et al., 2018). We use LSMs from the TRENDYv6 simulations for 1860–2015. To identify the contributions of CO₂ fertilization, climate, and land-use and land-cover change (LULCC) and management to the observed changes in SCA_NBP, we use outputs from three factorial simulations.

The models in simulation S3 were forced by (i) atmospheric CO₂ concentrations from ice core data and observations, (ii) historical climate reanalysis from the CRU-NCEP v8 (Viovy, 2016; Harris et al., 2014), and (ii) human-induced land-cover changes and management from a recent update of the land-use harmonization strategy (Hurtt et al., 2011) prepared for the next set of historical CMIP6 simulations, LUH2v2h (described below). Most models, though, still do not represent many of the management processes included in LUH2v2h. As summarized in Table A1 in Le Quéré et al. (2018), four models do not simulate wood harvest and three do not simulate cropland harvest. Two models simulate crop fertilization, tillage and grazing.

The models in simulation S2 were forced by (i) and (ii) described above with a fixed land-cover map from 1860. Simulation S2 estimates “natural” fluxes, and the difference between S2 and S3 outputs corresponds to anthropogenic CO₂ fluxes from LULCC. The models in simulation S1 were forced by changing atmospheric CO₂ and no climate change (recycling 1901–1920 values to simulate interannual variability) or LULCC. S1 thus provided changes in the terrestrial sink due to CO₂ fertilization, and the difference between S1 and S2 indicates the influence of climate change only. However, management practices (e.g. wood harvest), when simulated, are already included in S1 and S2 for some models. A baseline simulation with none of these effects (S0) was also performed to check for residual variability and trends. We selected only models providing spatially explicit outputs for the four simulations (S0, S1, S2 and S3) at monthly intervals (to evaluate seasonality; Table S1 in the Supplement).

We used NBP outputs selected for the period common to the inversion data, i.e. 1980–2015. NBP corresponds to the simulated net atmosphere–land flux (positive sign for a CO₂ sink), i.e. gross primary productivity (GPP) minus total ecosystem respiration (TER), fire emissions, and fluxes from LULCC and management (e.g. deforestation, agricultural and wood harvest, and shifting cultivation). All model outputs were resampled to a common regular latitude–longitude grid of 1^∘×1^∘.

2.3 Land cover and management

2.3.1 LUH2v2h

The LUH2v2h (Hurtt et al., 2011; available at http://luh.umd.edu/, last access: 7 July 2019) provides historical states and transitions of land use and management in a regular latitude–longitude grid of 0.25^∘×0.25^∘, covering the period 850–2015 at annual time intervals. Land-use states distinguish between primary and secondary natural vegetation (and forest and non-forest subtypes), managed pastures and rangelands, and multiple crop functional types. The updated data set includes several new layers of agricultural management, such as irrigation, nitrogen fertilization and biofuel management, and spatially explicit information about wood harvest constrained by LANDSAT data. Each LSM, however, may not simulate all the processes introduced in LUH2v2h, so the S3 results from each simulation might not be directly comparable.

2.3.2 ESA-CCI land cover

Land-cover information in LUH2v2h is combined with partial information on land use (e.g. rangeland in LUH2v2h can be either grassland or shrubland with low grazing disturbance). We therefore compared this information to annual land-cover maps at a latitude–longitude resolution of 0.5^∘×0.5^∘ based on the 300 m satellite-based land-cover data sets from ESA-CCI land cover (LC) (https://www.esa-landcover-cci.org/?q=node/175, last access: 1 August 2017) for 1992–2015. Data are provided for different vegetation types but were aggregated here for four main land-cover classes: forest, shrubland, grassland and cropland. The average distribution of these classes (forest and shrubland aggregated for readability) is shown in Fig. 2a. LUH2v2h was used for the statistical analysis of inversion and the LSM drivers (because it was the data set used to force the models), ESA-CCI data were used for the analysis of satellite-based vegetation data sets, and results were additionally compared with LUH2v2h.

2.4 Satellite-based vegetation data sets

We further evaluated trends in the activity and growth of vegetation for the different land-cover classes using three satellite-based data sets: leaf-area index (LAI), net primary production (NPP) and above-ground biomass (AGB) stocks. The LAI data set was calculated from satellite imagery from Global Inventory Modeling and Mapping Studies (GIMMS LAI3g) described by Zhu et al. (2015) for 1982–2015. LAI data were provided in two time steps per month on a regular latitude–longitude grid of 1∕12^∘ (subsequently aggregated to 0.5^∘). Smith et al. (2016) used the MODIS NPP algorithm and data for LAI and the fraction of photosynthetically active radiation from GIMMS to produce a 30-year global NPP data set, provided at monthly timescales for 1982–2011 at a latitude–longitude resolution of 1^∘×1^∘. The data are available at the NTSG data portal (https://wkolby.org/data-code/, last access: 4 August 2017). AGB stocks can be derived from estimates of vegetation optical depth derived from passive-microwave satellite measurements. Liu et al. (2015) produced a 20-year data set of AGB stocks for 1993–2012 based on measurements from a series of passive-microwave sensors. The data set is provided at a latitude–longitude resolution of 0.25^∘ ×0.25^∘ in annual time intervals and is available at http://www.wenfo.org/wald/global-biomass/ (last access: 13 February 2018). We tracked changes in LAI, NPP and AGB stocks for different land-cover types over time by selecting periods of at least 20 years common to ESA-CCI LC and the vegetation data sets (1992–2012 for LAI, 1992–2011 for NPP and 1993–2012 for AGB stocks). Vegetation variables were then aggregated for the four land-cover types at each time interval to account for land-cover changes.

3 Methods

3.1 Trends in seasonal-cycle amplitude (SCA_NBP)

The seasonal amplitude of CO₂ concentration is modulated by higher ecosystem CO₂ uptake during the growing season and increased emissions during the release period (TER) and thus controlled by the seasonal amplitude of NBP. We calculated SCA_NBP as the difference between peak uptake and trough for each year at the pixel scale, shown in Fig. 1a. However, since inversion fluxes have large uncertainty at the pixel level, we focused our analysis on SCA_NBP trends estimated from aggregated NBP over latitudinal bands or TransCom3 regions (Baker et al., 2006). Because we do not impose the timing of peak and trough, changes in SCA_NBP can be affected by the relative phase changes of GPP versus TER.

Figure 1 Variability in seasonal-cycle amplitude and trends from the inversions and LSMs. (a) Geographical distribution of SCA_NBP trends from the inversions (CAMS and CarboScope) and the multi-model ensemble mean (MMEM) from TRENDYv6 simulation S3 (all forcings). Both inversions estimated predominantly positive trends in SCA_NBP>40^∘ N (Fig. S1), so we defined two latitudinal bands, L_>40 N and L_25−40 N, for flux aggregation. (b, c) Aggregated SCA_NBP time series estimated by the inversions (CAMS in black and CarboScope s76 in grey) and S3 MMEM (red) for L_>40 N and L_25−40 N, respectively. The dashed lines indicate the linear fits used to calculate the slopes of the trends (corresponding colours), and the slopes and confidence intervals (95 %) are provided.

Figure 2 Regional distribution of the dominant land-cover types and SCA_NBP trends. (a) Land-cover map averaged over the study period for the three main land-cover classes (forest–shrubland, grassland and cropland) based on ESA-CCI annual land-cover data (1992–2015 average). (b) The continental regions correspond to the regions defined by Baker et al. (2006) and are delimited by bold lines: boreal and temperate North America (BorNA and TempNA), Europe (Eur), and boreal and temperate Eurasia (BorEA and TempEA). (c) Comparison of the SCA_NBP trends from the inversions to the trends estimated by the LSM experiments: S3, S2 (no LULCC) and S1 (no LULCC and no climate change). The bars for the inversions and LSMs indicate the average trend over each latitudinal band. The error bars for the inversions indicate the 95 % confidence levels for the trend values, and the vertical lines for the LSMs indicate interquartile ranges of the multi-model ensemble. The 95 % confidence interval for the MMEM was also calculated (see Methods).

The trend in SCA_NBP was calculated by a least-squares linear fit of annual values for 1980–2015, and confidence intervals were calculated based on Student's t distribution. We tested the robustness of estimated trends of inversions and LSMs for shorter periods by removing the first and last 1–10 years and trends of interannual variability by randomly removing 5 and 10 years of data 104 times. The significance of these trends was calculated using a Mann–Kendall test. We also compared different versions of CarboScope to evaluate the influence of the assimilated network size on the SCA_NBP trends (Fig. S1). We further calculated the trends for each of the sensitivity tests from CarboScope s85.

3.2 Process attribution

The three TRENDY experiments allow evaluating separately the effects of CO₂ fertilization, climate change, and LULCC in the models. The differences between S1 and S2 and between S2 and S3, however, could not isolate specific processes that may have contributed to the trend (e.g. cropland expansion versus afforestation or precipitation versus temperature). Furthermore, the LSMs may miss or simulate certain processes poorly that could influence SCA_NBP. Therefore, the attribution of drivers by the models is uncertain and should be cross-evaluated. Because inversions do not allow such partitioning between processes, a possible solution is to compare statistical attribution to drivers in inversions and LSMs.

We therefore compared the sensitivity of SCA_NBP estimated by the inversions and the LSMs by fitting a general linear model (GLM) using the iteratively reweighted least-squares method to eliminate the influence of outliers (Gill, 2000; Green, 1984). We tested the following variables (after unity-based normalization) as predictors: fertilization, irrigation, wood harvest, growing-season precipitation, growing-season temperature, atmospheric CO₂ concentration, and change in the extent of cropland and forest. These variables were taken from the corresponding data sets used to force TRENDYv6 models. All possible combinations of n predictors ($n=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{7}$) were tested, and for each value of n, the “best” model (according to Akaike’s information criterion) was chosen separately for each data set. Above n=4 no model showed improved fit compared to the models with fewer predictors. The coefficients from the GLM fit for each data set are shown in Fig. S5.

We further tested the robustness of the statistical relationships by fitting the GLM to the differences between each TRENDYv6 experiment. The significant predictors in the GLM fit to the LSMs in S3 should be detected in the corresponding factorial simulations, e.g. predictors associated with climate should be consistent for the fluxes estimated by the difference between S2 and S1 (effects of climate). The GLM fit to the partial fluxes for the effects of LULCC (S3−S2), climate (S2−S1) and CO₂ fertilization (S1−S0) are shown in Fig. S6.

4 Results

4.1 Large-scale patterns

4.1.1 Top-down estimates

Both inversions estimate increasingly positive trends in SCA_NBP with increasing latitude, even though CAMS shows heterogeneous patterns in North America, with strong decreasing trends for mid-latitudes (Figs. 1a and S1). Both inversions agree on significant positive SCA_NBP trends north of 40^∘ N (defined here as band L_>40 N) and non-significant trends for 25–40^∘ N (band L_25−40 N; Fig. 1b). In the L_>40 N band, CAMS and CarboScope s76 v4.1 estimate an SCA_NBP increase of 17.3±4.5 and 13.3±3.3 Tg C yr⁻², respectively. The uncertainties given for SCA_NBP trends represent here the uncertainty of the linear fit due to the year-to-year SCA_NBP variability (Methods). The difference between the CAMS and CarboScope inversions reflects part of the uncertainty in inversions due to their different choices in the ATM (including different atmospheric forcing and spatial resolution), the set of assimilated CO₂ data, the prior fluxes, and the a priori spatial and temporal correlation scales and is comparable to the uncertainty of the linear fit due to interannual variability.

This finding is corroborated by two further analyses of inversion uncertainties:

1. While both inversions assimilate atmospheric CO₂ measurements from Point Barrow, CAMS increasingly assimilates many other sites in the NH as they become available, helping to better constrain the CO₂ fluxes in mid-latitudes to high latitudes with time. Assimilating a non-stationary network of stations, however, possibly leads to spurious additional trends in SCA_NBP. To test this, we use different runs provided by CarboScope using more sites (but still fixed in number for each run; Fig. S2) for more recent periods. The results from CarboScope version s85 v4.1 (1985–2015) are generally consistent with CAMS, but version s93 v4.1 (1993–2015) estimates much stronger SCA_NBP trends (Table S1). A higher SCA_NBP trend in the period 1993–2015 is reported by both CAMS and CarboScope, which estimate very similar trends in L_>40 N (19.5 and 19.2 Tg C yr⁻², respectively).

2. CarboScope provides a set of sensitivity runs for s85 v4.1, varying some of the inversion's parameters (Fig. S3). Changes in the meteorological fields driving the transport model and the prior ocean fluxes have the largest effect on the SCA_NBP trends, giving L_>40 N trends of 8.6±4.9 Tg C yr⁻² (ERA-Interim instead of NCEP) and 13.9±5.6 Tg C yr⁻² (fixed ocean), respectively, both well within the uncertainty range (interannual variability affecting linear fit to SCA_NBP trend) estimated by the standard CarboScope s85 v4.1 (11.7±5.0 Tg C yr⁻²).

In summary, the ability of inversions to quantify the SCA_NBP trend is mostly limited by the intrinsic year-to-year SCA_NBP variability and less so by the amount of information available through the atmospheric data or by inversion settings.

4.1.2 Bottom-up estimates

The large-scale patterns of SCA_NBP trends from the LSM multi-model ensemble mean (MMEM) of simulation S3 (all forcings) are consistent with CarboScope inversion (Fig. 1a). The MMEM estimates are within the range of the inversions for most latitudes (Fig. S1) but always at the lower end of SCA_NBP trends reported by inversions. Consistent with inversions, LSMs report a significant trend in L_>40 N and a very weak (non-significant) trend in SCA_NBP in L_25−40 N (Fig. 1b). The overall MMEM trend in L_>40 N is significantly lower than in inversions (9.5±3.4 Tg C yr⁻², i.e. 55 %–71 % of inversions' estimates. The agreement between LSMs and inversions also varies depending on the period and set of inversions considered (LSMs capture 65 %–91 % of inversion trends in 1985–2015 and 74 %–75 % in 1993–2015; Table S1). The MMEM estimate for 1985–2015 (10.6±4.5 Tg C yr⁻²) is, in fact, even higher than the CarboScope inversion with different meteorological fields (8.6±4.9 Tg C yr⁻²). These results indicate that, despite a general underestimation of SCA_NBP trend in L_>40 N during 1980–2015 when compared to top-down estimates, the LSMs simulate the main spatiotemporal patterns in SCA_NBP trends consistent with inversions estimates, especially when accounting for the uncertainty in the latter.

To understand if recent improvements to the set of LSMs and their forcing in TRENDYv6 may have improved their performance in reproducing the SCA_NBP trend, we compared SCA_NBP trends from the previous intercomparison round (TRENDYv4). The MMEM from v6 estimates an SCA_NBP trend in L_>40 N that is 43 % higher than in than v4 (MMEM shown in Table S1 but evaluated for individual models). The specific reasons for improvement are hard to identify because of multiple model-dependent changes in the forcing, process simulation and parameterizations from v4 to v6 (Table 4 in Le Quéré et al., 2018).

In summary, we showed that the TRENDYv6 ensemble mean SCA_NBP trend captures the positive trends in the high latitudes and the lack of trend in the mid-latitudes given by inversions, and underestimates the magnitude of the high latitudes SCA_NBP trends by 9 %–45 %, depending on the inversion considered and period analysed.

4.2 Regional attribution

The comparison of SCA_NBP trends in large latitudinal bands may be useful in diagnosing general patterns but is less useful in diagnosing drivers of trends (e.g. climate and agriculture), since ecosystem composition, land management and climate effects are not necessarily separated along a latitudinal gradient. However, the comparison of inversions and models at the pixel scale is also not advisable because the sparse atmospheric network does not allow constraining the fluxes at this scale. We thus compared inversions and LSMs for the SCA_NBP trends over five sub-continental scale regions: boreal and temperate Eurasian and North American regions and Europe (“TransCom3” regions; Fig. 2). We then use LSMs for attributing SCA_NBP trends to different drivers using their factorial simulations (Methods).

Inversions and LSMs consistently attribute the increase in SCA_NBP mainly to boreal Eurasia, both in area-specific (Fig. 1a) and integrated values (Fig. 2b, 5.3–7.1 Tg C yr⁻² for inversions and 4.6 Tg C yr⁻² for MMEM, respectively), and to Europe (1.9–3.7 and 2.3 Tg C yr⁻²). The LSMs ascribed the trends in boreal Eurasia approximately equally to climate change and CO₂ fertilization (S1 and S2), with LULCC having a slight negative (i.e. decreasing) effect (compare S2 and S3), consistent with the results by Piao et al. (2017). In Europe, LSMs indicate negative contributions from both climate and LULCC. The negative effect of climate may be linked to increasingly drier conditions in this region (Greve et al., 2014) and to strong heatwaves in Europe in the early 21th century (Seneviratne et al., 2012). The negative contribution of LULCC indicated by LSMs in Europe does not support the idea that agricultural intensification or expansion drove an increase in SCA_NBP and is discussed further on. In temperate Eurasia, inversions disagree on the sign of SCA_NBP trends, and LSMs indicate weak positive trends dominated by the CO₂-fertilization effect. In boreal North America, LSMs estimate SCA_NBP trends very close to CarboScope estimates, mainly attributed to CO₂ followed by climate, whereas CAMS points to a trend close to zero because of cancelling regional trends with opposing sign (Fig. 1a). CAMS and CarboScope point to increasing SCA_NBP in temperate North America (1.4–1.6 Tg C yr⁻²), but the LSMs do not indicate any significant change (simulation S3). CAMS (which uses prior information with smaller a priori uncertainties than CarboScope together with a denser network) shows sharper regional differences than CarboScope, which illustrates that there are still substantial differences in the inversion at the scale of continental regions regarding SCA_NBP trends.

Aggregated over the two latitudinal bands (Fig. 2c), the MMEM indicates a dominant positive effect (increasing SCA_NBP) of CO₂ fertilization both in L_25−40 N and L_>40 N. In L_25−40 N, the CO₂ effect is offset by other factors: S1 differs significantly from S2 and S3, which have lower trends of SCA_NBP. In L_>40 N, the MMEM points to a positive effect of climate change in SCA_NBP trends, thus adding to the CO₂ effect. The MMEM suggest a negligible contribution of LULCC to the SCA_NBP trend in both latitudinal bands. The relative contributions of LULCC, climate and CO₂, however, differ between LSMs (Fig. S4). Most models nevertheless agree on non-significant SCA_NBP trends in L_25−40 N as well as on the predominant role of CO₂ fertilization and a non-significant contribution of LULCC to the trends in SCA_NBP in L_>40 N. Interestingly, models including carbon–nitrogen interactions had the weakest SCA_NBP trends (CABLE, ISAM and LPX-Bern), with the exception of CLM4.5, but we cannot draw conclusions from a small set of carbon–nitrogen models.

Figure 3 Statistical attribution of drivers of SCA_NBP estimated by the inversions and LSMs. The main drivers of SCA_NBP are presented for (a) L_>40 N and (b) L_25−40 N and are calculated as the product of the coefficients of a general linear model fit on SCA_NBP using a number of predictors (normalized) and their corresponding trends. Fertilization, irrigation, wood harvest, growing-season precipitation, growing-season temperature and atmospheric CO₂ concentration were tested as predictors, and the best fit was chosen for each data set: CAMS (dark grey), CarboScope s76 (light grey) and the MMEM (red). The bars indicate the contribution of each predictor to the trend in SCA_NBP, error bars indicate the corresponding 95 % confidence intervals and the symbols indicate significant GLM fits (two asterisks, one asterisk and crosses indicate p<0.01, p<0.05 and p<0.1, respectively).

4.3 Driving processes

Figure 3 shows the relative contributions of the predictors (weighted by their trends) found to SCA_NBP trends in both latitudinal bands. The coefficients of the GLM fit are shown in Fig. S5.

The GLMs provide a better fit the trend of SCA_NBP in L_>40 N (explaining 57 %–74 % of the variance; Table S2) than for L_25−40 N (8 %–49 % only). The GLM fit to inversions and to the MMEM identified CO₂ fertilization as the most important factor explaining (statistically) the SCA_NBP trends in both latitudinal bands, consistent with S1 (Figs. 1, S4, S5 and S6), although the CO₂-fertilization effect was weaker for the GLM fit to LSMs than for inversions in region L_>40 N. The statistical models for inversions and LSMs agreed on a significant negative contribution of warming in both latitudinal bands but a stronger contribution in L_25−40 N. GLM models fitted to LSMs and CarboScope also point to changes in forest area contributing to increased SCA_NBP, and changes in crop area have a negative effect in SCA_NBP from LSMs. In L_25−40 N, the GLM fit to LSMs further points to a small negative contribution of wood harvest to SCA_NBP trends, and, for CAMS, negative effects of irrigation and fertilization are also significant. The statistical attribution of SCA_NBP trends in LSMs is generally consistent with the factorial simulations but is mostly clearer for the CO₂-fertilization effect than for the other drivers. In the difference between factorial simulations (Fig. S6), some drivers appear to have strong interactive effects, e.g. the effect of CO₂ is significantly negative for S3−S2 (LULCC). This could be explained by higher emissions from LULCC under higher CO₂ concentrations from the loss of additional sink capacity (Pongratz et al., 2014). The key role of CO₂ fertilization in the observed changes is in line with Piao et al. (2017), but our results challenge some of the previously proposed hypotheses to account for the increase in seasonal CO₂ exchange, as addressed below.

5 Discussion

5.1 Confronting hypotheses

5.1.1 Contribution of LULCC

Agricultural intensification and expansion occurred mainly in latitudes below 45^∘ N (Gray et al., 2014), and inversions and LSMs reported instead a peak in the amplitude of land-surface CO₂ exchange for latitudes above 45^∘ N (Figs. 1 and S1). Furthermore, our regional attribution identifies Eurasia as the region contributing most to increasing SCA_NBP; this region is dominated by natural ecosystems (Fig. S3) and has experienced very little land use change (Verburg et al., 2015) over the previous decades. Additionally, factorial LSM simulations indicate a negligible contribution of LULCC and management to SCA_NBP trends not only at the latitudinal-band scale but also regionally (Figs. 1c and 2).

This, though, could not in itself falsify the hypothesis that agricultural intensification is a key driver of SCA_NBP trends because most LSMs still do not include processes that could intensify cropland NPP over time, such as better cultivars, fertilization and irrigation. Still, management practices are neither a significant predictor for GLM fitted to LSMs nor for inversions, with the exception of CAMS. CarboScope further identifies a negative effect of cropland expansion to SCA_NBP in L_>40 N rather than a positive one, which partly challenges the contribution of cropland expansion (Gray et al., 2014) to SCA_NBP. Our results are consistent with those by Smith et al. (2014) that show that NPP generally decreased following conversion from natural ecosystems to cropland, except in areas of highly intensive agriculture, such as the midwestern USA. Increasing crop productivity (intensification) could partly explain increasing SCA_NBP. However, satellite-based data for LAI (Zhu et al., 2016), NPP (Smith et al., 2016) and above-ground biomass (AGB) carbon stocks (Liu et al., 2015) for different land-cover classes from ESA-CCI LC (Fig. S7) indicate that the increase in crop productivity accounted for only a small fraction of the hemispheric trends in ecosystem productivity, consistent with the crop productivity stagnation in Europe and Asia identified by Grassini et al. (2013). We also compared the trends in remote-sensing variables for land-cover classes from LUH2v2h, with similar results.

Previous studies suggesting a large role of the green revolution in SCA trends focused on a longer period, starting in the 1960s. The acceleration of SCA_NBP reported by inversions and LSMs (Table S1) concurrent with crop productivity stagnation indicates that since the 1980s agriculture intensification is not likely to be the main driver of the increase in SCA. Even in the intensive agricultural areas in the midwestern USA, CAMS estimates contrasting negative and positive trends (Figs. 1a, S8). Eddy-covariance flux measurements (only for 7–13 years) in the areas of intensive agriculture in the USA show a weak relationship between trends in NBP and trends in SCA_NBP, showing mostly non-significant trends in SCA_NBP (Fig. S8).

5.1.2 Contribution of warming

We found that warming during the growing season had a negative effect on SCA_NBP trends in both latitudinal bands, although this effect is uncertain for LSMs in L_>40 N. Annual temperature used in the statistical models was also negatively correlated with SCA_NBP, but the correlation was only significant for CAMS.

The negative relationship with growing-season temperature (T) at the mid-latitudes may be explained by warmer temperature increasing atmospheric demand for water (Novick et al., 2016) and inducing soil-moisture deficits in water-limited regions in summer (Seneviratne et al., 2010) or increased fire risk (Peñuelas et al., 2017) that reduces the summer minimum of SCA_NBP. This negative effect of temperature can explain the negative contribution of climate to the simulated SCA_NBP trends in L_25−40 N given by the factorial simulations (Fig. 2c).

The negative statistical relationship found between the trend of SCA_NBP and T in L_>40 N challenges the assumption that a warming-related increase in plant productivity in high latitudes necessarily increases the seasonal CO₂ exchange (Keeling et al., 1996; Graven et al., 2013; Forkel et al., 2016). Although the MMEM shows a small positive contribution of climate in L_>40 N, LSMs diverge on the contribution of climate in this latitudinal band (Fig. S4). Moreover, the factorial simulations in Fig. 2c allow evaluating the impact of changes in all climate variables (e.g. also rainfall and radiation) in addition to temperature. A negative relationship between T and SCA has, though, also been reported by Schneising et al. (2014) for interannual changes in the SCA of total column CO₂ for 2004–2010. Yin et al. (2018) have further shown that, at latitudes between 60^∘ N and 80^∘ N, the relationship between SCA and T has transitioned from positive in the early 1980s to negative in recent decades, reconciling the results by Keeling et al. (1996) and Schneising et al. (2014).

In Fig. S9 we present a conceptual scheme of the impacts of warming in SCA_NBP through its component fluxes. Generally, warming in high latitudes has been associated with a longer growing season and increased GPP (Piao et al., 2008), which would contribute to increased SCA_NBP through increased productivity during the “uptake period” and increased decomposition (due to more litter) during the “release period”. However, a weakening of this relationship has been reported (Piao et al., 2014; Peñuelas et al., 2017). Other processes can, though, contribute to the negative relationship between SCA_NBP and T reported here and in other studies. The empirical negative relationship between trends in SCA_NBP and warming at the higher latitudes may be due to (i) a stronger effect of T on total ecosystem respiration (TER) than on GPP during the uptake period, (ii) a negative response of ecosystem productivity to warming during the uptake period or (iii) indirect negative effects of T on decomposition during the release period.

Figure 4 Emerging relationships between LSM sensitivities to climate and CO₂ and their SCA_NBP trends. The SCA_NBP trend for L_>40 N estimated by each inversion (grey intervals) and corresponding responses of SCA_NBP to (a) T and (b) CO₂ (as calculated in Fig. 3 but considering the scores of the regression only; shown in Fig. S5) are compared to the results from individual models (simulation S3; coloured markers). In (c) the SCA_NBP trend for L_>40 N estimated by individual models is compared with the simulated sensitivity of TER to T. The shaded areas indicate the inversion ranges, and the distribution of the grey lines shows uncertainty in the relationship between each pair of variables.

Evidence nevertheless supports negative effects of warming on SCA trends. Temperature increase in recent decades has been associated with widespread reduction in the extent and depth of snow cover (Kunkel et al., 2016) and in the number of days with snow cover (Callaghan et al., 2011). Snow has an insulating effect, so snow-covered soil during winter can be kept at relatively constant temperatures, several degrees above the air temperature (>10 ^∘C), which promotes respiration of soil C (Nobrega and Grogan, 2007). Soils become subject to more fluctuations in temperature, and become colder, as the snow cover recedes or becomes thinner. Yu et al. (2016) reported that respiration suppression due to a reduction in snow cover in winter may account for as much as 25 % of the increase in the annual CO₂ sink of northern forests. A decrease in respiration in response to warming during the release period could thus decrease SCA_NBP, but the effect of growing-season temperature was stronger in our study. The expansion of vegetation in Arctic tundra, particularly shrubland, has been linked to warming trends but also depends on soil-moisture and permafrost conditions (Elmendorf et al., 2012). Many regions of dry tundra and low arctic shrubland (Walker et al., 2005) experience summer drought or soil-moisture limitations even though northern regions are usually considered to be energy-limited (Greve et al., 2014). Indeed, Myers-Smith et al. (2015) found a strong soil-moisture limitation of the (positive) sensitivity of shrub growth to temperature in summer, possibly associated with the limitation of growth due to drought and/or with reduced growth and dieback due to standing water during thawing. CAMS indicates a decrease in SCA_NBP in eastern regions in boreal North America (Fig. 1a), where Myers-Smith et al. (2015) reported negative sensitivity of shrub-growth to temperature. The coarse network and large correlation lengths used by CarboScope do not allow such regional contrasts to be resolved. Most process-based models lack a detailed representation of processes described above – e.g. a realistic effect of snow insulation on soil temperatures, soil freezing and thawing (Koven et al., 2009; Peng et al., 2016; Guimberteau et al., 2017) – potentially overestimating the net sink response to temperature changes (Myers-Smith et al., 2015). Moreover, soil-moisture limitation due to temperature increase could also contribute to decreased TER by limiting microbial activity, which is currently not simulated in most LSMs. This may in turn explain why LSMs underestimate the negative effect of temperature in SCA_NBP in the high latitudes compared to CAMS (Figs. 3 and S4).

5.2 Evaluating model biases

Wenzel et al. (2016) proposed that the observed sensitivity of SCA_NBP to CO₂ was an emergent constraint on future terrestrial photosynthesis, but their study focused on simulations by an earth-system model that excluded the effects of climate change (i.e. the radiative feedback of CO₂ to climate was not considered). Our results are consistent with a strong increase in the peak uptake due to the effect of CO₂ fertilization driven by GPP as proposed by Wenzel et al. (2016). The negative effect of temperature in our study (Fig. 3), although weaker than the positive effect of absolute CO₂ concentration, suggested that warming partly cancelled out the increase in SCA_NBP expected from the effect of fertilization alone. We propose that other processes partly control SCA_NBP trends linked to reduced decomposition under lower snow cover (Yu et al., 2016) or to emerging limitations to growth in response to water limitation (Elmendorf et al., 2012; Myers-Smith et al., 2015). Additionally, while the sensitivity of productivity to the CO₂-fertilization effect is expected to decrease, the control of respiration by temperature should increase non-linearly (Piao et al., 2014, 2017; Peñuelas et al., 2017), suggesting a progressively dominant (negative) influence of warming on SCA_NBP. The degree of such an offset would likely depend on the thresholds of soil temperature and water limitation that are complex and thus difficult to assess and require process-based modelling. Our results imply that future constraints of productivity based only on the CO₂ effect (as in Wenzel et al., 2016) may overestimate future GPP.

We evaluated whether the differences between the observed SCA_NBP trends (significant only in L_>40 N) and those simulated by the LSMs could be associated with the modelled sensitivities to atmospheric CO₂ concentration (CO₂) and growing-season temperature (T) in L_>40 N (Fig. 4). In Fig. 4, only models with too small a sensitivity of SCA_NBP to T produce a realistic trend of SCA_NBP. In contrast, the models indicating sensitivities to T and CO₂ more similar to those estimated by the inversions tend to underestimate the trend in SCA_NBP.

Why are the LSM sensitivities of SCA_NBP to T positively correlated with their long-term SCA_NBP trend (Fig. 4) even though CO₂ is a stronger driver of the simulated SCA_NBP trend (Fig. 3)? We found a clear relationship between the model bias in the trend of SCA_NBP and the sensitivity to CO₂ fertilization in S3 (in line with Wenzel et al., 2016), but we also found a compensatory effect, where models that overestimate the sensitivity of SCA_NBP to T tend to underestimate the sensitivity to CO₂ and vice versa. LSMs tend to overestimate the sensitivity of SCA_NBP to T and underestimate the sensitivity to CO₂, compared to the observation-based constraints from inversions. LSMs often compensate too strong (or too weak) a simulated water stress or temperature sensitivity by adjusting photosynthesis parameters (that control CO₂ fertilization) during model optimization to match the observed net terrestrial sink. This compensatory effect has previously been reported by Huntzinger et al. (2012) for the mean terrestrial sink; we find that it could also affect the trends in seasonal CO₂ exchange.

We argue that the trend of SCA_NBP can differ between models due to (a) differences in their NPP response to T and CO₂, (b) differences in turnover times of short-lived C pools by which increased NPP is coupled to increased winter decomposition, and (c) phase shifts between GPP and ecosystem respiration. These phase shifts may be associated with errors in the phase and amplitude of simulated ecosystem respiration, arising from factors such as (i) representing soil carbon stocks as pools with discrete turnover times and associated effective soil depths (Koven et al., 2009) and (ii) neglect of seasonal acclimation effects on autotrophic and heterotrophic respiration. The sensitivities of NPP to CO₂ and T between models are strongly and consistently correlated with the compensatory effect of the model parameterizations (Fig. S10), but we find no clear relationship between the biases of the modelled SCA_NBP trend and the sensitivity of NPP to T (Fig. S11), suggesting a key role of respiration. Indeed, the models with SCA_NBP trends closer to observations tend to be associated with a lower sensitivity of ecosystem respiration to growing-season temperature (Fig. S4c). Too large a turnover of short-lived pools in a model should produce too small an increase in the SCA_NBP amplitude (i.e. increased respiration during the uptake period followed by too little during the release period) for a given sensitivity of NPP to CO₂ or climate. A recent study by Jeong et al. (2018) has reported that ecosystem carbon-cycle models (not used in this study) underestimated changes in carbon residence times in northern Alaska. The evaluation of the effect of model turnover times in SCA_NBP requires a deeper analysis of transfers between litter and soil organic carbon pools and can be verifiable in future simulations.

6 Conclusions

Based on our assessment of atmospheric observations and the most advanced land-surface model simulations, the most likely explanation of the seasonal cycle of atmospheric CO₂ at high latitudes is the CO₂ fertilization of photosynthesis in unmanaged high-latitude ecosystems, especially in the Eurasian boreal forests. Our study further points to key processes that need to be developed to better simulate NBP responses to changing climate, especially to Arctic warming, in particular productivity limitations and the decomposition terms. Our results indicated that the signal of the SCA_NBP trend contains valuable information for the turnover times of short-term pools, which await further investigation.