Contrasting effects of CO2 fertilisation, land-use change and warming on seasonal amplitude of northern hemisphere CO2 exchange

Continuous atmospheric CO2 monitoring data indicate an increase in seasonal-cycle amplitude (SCA) of CO2 exchange in northern high latitudes. The major drivers of enhanced SCA remain unclear and intensely debated with land-use change, CO2 fertilization and warming identified as likely contributors. We integrated CO2-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 of CO2-fluxes for 1980–2015. The LSMs generally 5 reproduce the latitudinal increase in SCA trends within the inversions range. Inversions and LSMs attribute SCA increase to


5
The increase in the amplitude of seasonal atmospheric CO 2 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 2 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  and appears to be driven primarily by changes in seasonal growth dynamics of terrestrial ecosystems (i.e., net 10 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 2 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 15 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 2 concentrations using atmospheric transport models (ATM) 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 20 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 2 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 in other sites, LSMs generally captured the sign of the 25 trend in SCA but either under-or over-estimated its magnitude. Piao et al. (2017) further reported that CO 2 fertilisation and climate change drove the increase in SCA for sites >50°N, but that at mid-latitude sites land use, oceanic fluxes, fossil-fuel emissions, as well as trends in atmospheric transport may have contributed to the SCA trends.
Atmospheric inversions provide a consistent framework for assimilating in-situ CO 2 concentration observations to estimate net CO 2 surface fluxes while accounting for errors in the prior fluxes and for some errors in the ATM (Peylin et al., 2013). SCA is fully consistent with the amplitude of CO 2 concentrations in all stations of the observational network used and constitute a direct benchmark for SCA N BP 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 2 fluxes at northern extra-tropical latitudes between 1980-2015 to: (i) assess the ability of those LSMs to simulate inversion-based trends in SCA N BP ; 5 (ii) attribute the trends in SCA N BP to specific regions in the Northern Hemisphere; (iii) attribute the relative importance of drivers using the ensemble model framework.
The two inversions used here -the Copernicus Atmosphere Monitoring Service (CAMS) inversion system (Chevallier et al., 2005) and the Jena CarboScope inversion (Rödenbeck, 2005) -solve for fluxes on their ATM grid, thus minimising aggregation errors for large regions (Kaminski and Heimann, 2001). The two inversions differ in a number of characteristics (see Methods), 10 particularly the ATM, the prior information and the observations assimilated: CAMS includes a time-varying number of multi- year air-sampling sites as they are available (constraining better spatial patterns), and CarboScope keeps a fixed set of sites covering a given period (avoiding artefacts in the time series related to the appearance or disappearance of measurement sites).
We use NBP simulated by a set of 11 LSMs from the recent TRENDYv6 inter-comparison (Le Quéré et al., 2018) in three distinct experiments where models are forced with: changing CO 2 only (S1); CO 2 and changing climate (S2); changing 15 climate, CO 2 , and LULCC (S3). Some of the TRENDYv6 models used here simulate some key management processes (4 models include wood harvest, 1 model irrigation, and 2 fertilization) and more models now include nitrogen cycling, and improved soil processes, which were missing in previous intercomparisons (Arneth et al., 2017)  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 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 inversions systems that provide results for the study period between 1980 and 2015. We analysed 25 monthly surface CO 2 fluxes estimated by the inversion systems from the Copernicus Atmosphere Monitoring Service (CAMS) (Chevallier et al., 2005(Chevallier et al., , 2010 and from Jena CarboScope (Rödenbeck et al., 2003;Rödenbeck, 2005).
Here we use CAMS version r16v1 (http://atmosphere.copernicus.eu/) selected for the period between 1980-2015, that provides estimates of ocean and terrestrial fluxes at 1.9°latitude ×3.75°longitude resolution. The CAMS inversion system assimilates observations from a variable number of atmospheric CO 2 monitoring sites (119 in total providing at least 5 years 30 of measurements) and uses the transport model from the LMDz General Circulation Model (LMDz5A) nudged to ECMWFanalysed 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) provides several versions that assimilate a temporally consistent set of observations. We used these versions for the study period  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 5 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: "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; "short" 10 -reducing the length of a-priori spatial correlations. The resulting latitudinally-integrated SCA N BP and respective trends are shown in Figure S2.

Land-surface Models
Land-surface models (LSMs) provide a bottom-up approach to evaluate terrestrial CO 2 fluxes (i.e. net biome productivity, NBP), and allow deeper insight into the mechanisms driving changes in C-stocks and fluxes. The TRENDY intercomparison 15 project compiles simulations from state-of-the-art LSMs to evaluate terrestrial energy, water and CO 2 exchanges since the preindustrial period (Sitch et al., 2015;Le Quéré et al., 2018). We use LSMs from the TRENDY v6 simulations for 1860-2015. To identify the contributions of CO 2 fertilisation, climate, and LULCC and management to the observed changes in SCA N BP , we use outputs from three factorial simulations.
The models in simulation S3 were forced by (i) atmospheric CO 2 concentrations from ice core data and observations, (ii) 20 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 (Hurtt et al., 2011) prepared for the next set of historical CMIP6 simulations, LUH2v2h (described below). Most models still do not represent many of the management processes included in LUH2v2h, though. 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. 25 The models in simulation S2 were forced by (i) and (ii) with fixed land-cover map from 1860. Simulation S2 estimates "natural" fluxes, and the difference between S2 and S3 outputs corresponds to anthropogenic CO 2 fluxes from LULCC. The models in simulation S1 were forced by changing atmospheric CO 2 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 2 fertilisation, and the difference between S1 and S2 indicates the influence of climate change only. However, management practices (e.g. wood-30 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, Supplementary Table 1).
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 2 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°. such as irrigation, nitrogen fertilisation, 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.

ESA-CCI Land-Cover
Land-cover information in LUH2v2h is combined with partial information on land use (e.g. rangeland in LUH2v2h can be 15 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 LC (https://www.esa-landcover-cci.org/?q=node/175) for 1992-2015. Data are provided for different vegetation types, but here were aggregated 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 Figure 2a. LUH2 was used for the statistical analysis of 20 inversion and the LSM drivers (because it was the data set used to force the models), and ESA-CCI data were used for the analysis of satellite-based vegetation data sets.

Satellite-based vegetation datasets
We further evaluated trends in the activity and growth of vegetation for the different land-cover classes using three satellitebased data sets: leaf-area index (LAI), net primary production (NPP), and aboveground biomass (AGB) stocks. The LAI data 25 set was calculated from satellite imagery from Global Inventory Modeling and Mapping Studies (GIMMS LAI3g) described by  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 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°i n annual time intervals and is available at http://www.wenfo.org/wald/global-biomass/ (last access 13/02/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 (1992for LAI, 1992-2011for NPP, and 1993 for AGB stocks).

5
Vegetation variables were then aggregated for the four land-cover types at each time interval to account for land-cover changes.

Trends in seasonal-cycle amplitude (SCA)
The seasonal amplitude of CO 2 concentration is modulated by higher ecosystem CO 2 uptake during the growing season and increased emissions during the release period (TER) and thus controlled by the seasonal amplitude of NBP. We calculated 10 SCA N BP as the difference between peak uptake and trough for each year, at pixel scale shown in Figure 1a. However, since inversion fluxes have large uncertainty at pixel-level we focused our analysis on SCA 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 N BP can be affected by the relative phase changes of GPP versus TER.
The trend in SCA N BP was calculated by a least-squares linear fit of annual values for 1980-2015, and confidence intervals 15 were calculated based on the 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 N BP trends ( Figure S1). We further calculated the trends for each of the sensitivity tests from CarboScope s85.

Process attribution
The three TRENDY experiments allow evaluating separately the effects of CO 2 fertilisation, 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 poorly certain processes that could influence SCA N BP . Therefore, the attribution of drivers 25 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 N BP 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 normalisation) as predictors: fertilization, irrigation, wood harvest, forest. These variables were taken from the corresponding datasets used to force TRENDYv6 models. All possible combinations of n predictors (n=1,2, . . . , 7) were tested, and for each value of n, the "best" model (according to Akaike's information criterion) was chosen separately for each dataset. Above n=4 no model showed improved fit compared to the models with less predictors. The coefficients from the GLM fit for each dataset are shown in Figure S.
We further tested the robustness of the statistical relationships by fitting the GLM to the differences between each TRENDYv6 5 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 2 fertilisation (S1-S0) are shown in Figure S5. TgC.yr -2 and 13.3±3.3 TgC.yr -2 , respectively. The uncertainties given for SCA N BP trends represent here the uncertainty of the linear fit due to the year-to-year SCA N BP 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 2 data, the prior fluxes, and the a-priori spatial and temporal correlation 20 scales, and is comparable to the uncertainty of the linear fit due to inter-annual variability. This finding is corroborated by two further analyses of inversion uncertainties:

Results and discussion
(1) While both inversions assimilate atmospheric CO 2 measurements from Point Barrow, CAMS increasingly assimilates many other sites in the NH as they become available, helping to better constrain the CO 2 fluxes in mid-to high-latitudes with time. Assimilating a non-stationary network of stations, however, possibly leads to spurious additional trends in SCA N BP . 25 To test this, we use different runs provided by CarboScope using more sites (but still fixed in number for each run, Figure S2) for more recent periods. The results from CarboScope version s85 v4.1  are generally consistent with CAMS, but version s93 v4.1 (1993-2015) estimates much stronger SCA N BP trends (Table S1). A higher SCA N BP trend in the period 1993-2015 is reported by both CAMS and CarboScope, which estimate very similar trends in L >40N (19.5 TgC.yr -2 and 19.2 TgC.yr -2 respectively).
In summary, the ability of inversions to quantify the SCA N BP trend is mostly limited by the intrinsic year-to-year SCA N BP variability, less so by the amount of information available through the atmospheric data or by inversion settings. To understand if recent improvements to the set of LSMs and their forcing in TRENDYv6 may have improved their performance in reproducing the SCA N BP trend, we compared SCA N BP trends from the previous intercomparison round (TRENDYv4). The MMEM from v6 estimates an SCA N BP trend in L >40N 43% higher than in than v4 (MMEM shown in 20 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 N BP trend captures the positive trends in the high latitudes and the lack of trend in the mid-latitudes given by inversions, and under-estimates the magnitude of the high latitudes 25 SCA N BP trends by 9-45%, depending on the inversion considered and period analysed.

Regional attribution
The comparison of SCA N BP trends in large latitudinal bands may be useful in diagnosing general patterns, but is less useful to diagnose drivers of trends (e.g. climate, 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 pixel scale is also America regions, and Europe ("TransCom3" regions, Figure 2). We then use LSMs for attributing SCA N BP trends to different drivers using their factorial simulations (Methods).
Inversions and LSMs consistently attribute the increase in SCA N BP mainly to boreal Eurasia, both in area specific ( Figure   1a) and integrated values (Figure 2b, 5.3-7.1 TgC.yr -2 for inversions and 4.6 TgC.yr -2 for MMEM, respectively) and to Europe (1.9-3.7 TgC.yr -2 and 2.3 TgC.yr -2 ). The LSMs ascribed the trends in boreal Eurasia approximately equally to climate change 5 and CO 2 fertilisation (S1 and S2), with LULCC having a slight negative (i.e. decreasing) effect (compare S2 and S3), consistent with the results by . 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 21 th 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 N BP and 10 is discussed further on. In temperate Eurasia, inversions disagree on the sign of SCA N BP trends and LSMs indicate weak positive trends dominated by the CO 2 fertilisation effect. In boreal North America, LSMs estimate SCA N BP trends very close to CarboScope estimates, mainly attributed to climate, whereas CAMS points to a trend close to zero because of cancelling regional trends with opposing sign (Figure 1a). CAMS and CarboScope point to increasing SCA N BP in temperate North America (1.4-1.6 TgC.yr -2 ), but the LSMs do not indicate any significant change (simulation S3). CAMS (which uses prior 15 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 N BP trends.
Aggregated over the two latitudinal bands (Figure 2c), the MMEM indicates a dominant positive effect (increasing SCA N BP ) of CO 2 fertilization both in L 25−40N and L >40N . In L 25−40N , the CO 2 effect is offset by other factors: S1 differs significantly 20 from S2 and S3, which have lower trends of SCA N BP . In L >40N , the MMEM points to a positive effect of climate change in SCA N BP trends, thus additive to the CO 2 effect. The MMEM suggest a negligible contribution of LULCC to the SCA N BP trend in both latitudinal bands. The relative contributions of LULCC, climate and CO 2 however, differ between LSMs ( Figure   S4). Most models nevertheless agree on non-significant SCA N BP trends in L 25−40N as well as on the predominant role of CO 2 fertilisation and a non-significant contribution of LULCC to the trends in SCA N BP in L >40N . Interestingly, models 25 including carbon-nitrogen interactions had the weakest SCA N BP trends (CABLE, ISAM and LPX-Bern), excepting CLM4.5 but we cannot draw conclusions from a small sent of carbon-nitrogen models.

Driving processes
Validating the attribution of SCA N BP trends to CO 2 , climate, and LULCC by the LSMs at large scale is challenging. However, insights into the consistency of SCA N BP drivers can be obtained by statistical analysis. We fitted multiple general linear 30 regression models (GLM) to the SCA N BP from the inversions and the S3 MMEM over each latitudinal band using multiple predictors from the TRENDYv6 forcing datasets (CO 2 concentration, climate variables and changes in land-cover and management practices). For each dataset, we identified the statistical model that could best explain SCA N BP trends with least number of predictors (Table S2). We then tested the results from this statistical attribution for the MMEM against the corresponding  Figure 3 shows the relative contributions of the predictors (weighted by their trends) found to SCA N BP trends in both latitudinal bands. The coefficients of the GLM fit are shown in Figure S5.
The GLMs provide a better fit the trend of SCA N BP in L >40N (57-74% of the variance, Table S2) than for L 25−40N (8-49% only). The GLM fit to inversions and to the MMEM identified CO 2 fertilisation as the most important factor explaining (statistically) the SCA N BP trends in both latitudinal bands, consistent with S1 (Figures 1, S4, S5 and S6), even though the 5 CO 2 fertilization effect was weaker for the GLM fit to LSMs than for inversions in region L >40N . The statistical models for inversions and LSMs agreed on a significant negative contribution of warming in both latitudinal bands, but stronger in L 25−40N . In L >40N , GLM models fitted to LSMs and CarboScope also point to changes in forest area contributing to increase SCA N BP , and changes in crop area have a negative effect in SCA N BP from LSMs. In L 25−40N , the GLM fit to LSMs further points to a negative contribution of wood-harvest to SCA N BP trends, and the fit to CAMS a weak negative effects of irrigation 10 and fertilization. The statistical attribution of SCA N BP trends in LSMs is consistent with the factorial simulations, although the negative effect of temperature is only significant in L 25−40N . The key role of CO 2 fertilisation in the observed changes is in line with Piao et al. (2017), but our results challenge some previously proposed hypotheses to account for the increase in seasonal CO 2 exchange, as addressed below. LSMs reported instead a peak in the amplitude of land surface CO 2 exchange for latitudes above 45°N (Figure 1 and S1). Furthermore, our regional attribution identifies Eurasia as the region contributing most to increasing SCA; this region is dominated by natural ecosystems ( Figure S3) and has experienced very little land use change (Verburg et al., 2015) over the past decades. 20 Additionally, factorial LSM simulations indicate a negligible contribution of LULCC and management to SCA N BP trends at latitudinal-band scale but also regionally (Figures 1c and 2). This, though, could not in itself falsify the hypothesis that agricultural intensification is a key driver of SCA N BP trends, because most LSMs still do not include processes that could intensify cropland net primary productivity (NPP) over time such as better cultivars, fertilization, irrigation. Still, management practices are not a significant predictor for GLM fitted to 25 LSMs, but also not for inversions, excepting CAMS. CarboScope further identifies a negative effect of cropland expansion to SCA N BP in L >40N rather than a positive one, which partly challenges the contribution of cropland expansion (Gray et al., 2014) to SCA N BP . Our results are consistent with those by Smith et al. (2014) that show that net primary productivity (NPP) generally decreased following conversion from natural ecosystems to cropland, except in areas of highly intensive agriculture such as midwestern USA. Increasing crop productivity (intensification) could partly explain increasing SCA N BP . However, 30 satellite-based data for LAI (Zhu et al., 2016), NPP (Smith et al., 2016)  trends in ecosystem productivity, consistent with crop productivity stagnation in Europe and Asia identified by Grassini et al. (2013).

Confronting Hypoteses
Previous studies suggesting a large role of the green revolution in SCA N BP trends have focused on a longer period, starting in the 1960s. The acceleration of SCA N BP reported by inversions and LSMs (Table 1) 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.

5
Even in the intensive agricultural areas in the US Midwest, CAMS estimates contrasting negative/positive trends (Figure 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 N BP , showing mostly non-significant trends in SCA N BP ( Figure S8).

Contribution of warming
We found that warming during the growing season had a negative effect on SCA N BP trends in both latitudinal bands, although 10 this effect is uncertain for LSMs in L >40N . Annual temperature used in the statistical models was also negatively correlated with SCA N BP , 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  that reduce the summer minimum of SCA N BP . The negative statistical relationship found between 15 the trend of SCA N BP and T in L >40N challenges the assumption that warming-related increase in plant productivity in highlatitudes necessarily increases the seasonal CO 2 exchange (Keeling et al., 1996;Graven et al., 2013;Forkel et al., 2016). Such a negative relationship, however, has also been reported by (Schneising et al., 2014) for interannual changes in the SCA N BP of total column CO 2 for 2004-2010. Yin et al. (2018) have further shown that, at latitudes between 60°N and 80°N, the relationship between SCA NBP and T has transitioned from positive in the early 1980s, to negative in recent decades, reconciling the results 20 by Keeling et al. (1996) and Schneising et al. (2014).
The empirical negative relationship between trends in SCA N BP and warming at the higher latitudes may be due to either (i) indirect negative effects of T on decomposition during the "release period"; (ii) a negative response of ecosystem productivity to warming during the "uptake period"; (iii) a stronger effect of T on total ecosystem respiration (TER) than on GPP during the growing season. Some of these effects are counter-intuitive, because warming in high-latitudes is usually associated with 25 longer growing-season and increased GPP (Piao et al., 2008), although a weakening of this relationship has been reported (Piao et al., 2014;Peñuelas et al., 2017).
Evidence nevertheless supports negative effects of warming on SCA trends. Temperature increase in recent decades has been associated with widespread reduction in 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 30 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 2 sink of northern forests. A decrease in respiration in response to warming during the release period could thus decrease SCA N BP , 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) 5 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 N BP in eastern regions in boreal North America (Figure 1a), where Myers-  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 10 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, soilmoisture limitation due to temperature increase could also contribute to decrease 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 N BP in the high-latitudes compared to CAMS (Figure 3 and S4).  Wenzel et al. (2016) proposed that the observed sensitivity of SCA N BP to CO 2 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 2 to climate was not considered). Our results are consistent with a strong increase in the peak uptake due to the effect of CO 2 fertilisation driven by gross primary production (GPP) as proposed by Wenzel et al. (2016). The 20 negative effect of temperature in our study (Figure 3), although weaker than the positive effect of absolute CO 2 concentration, suggested that warming partly cancelled out the increase in SCA N BP expected from the effect of fertilisation alone. We propose that other processes partly control SCA N BP 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 2 fertilisation effect is expected to decrease, whereas the 25 control of respiration by temperature should increase nonlinearly (Piao et al., 2014Peñuelas et al., 2017), suggesting a progressively dominant (negative) influence of warming on SCA N BP . 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 2 effect (as in Wenzel et al. (2016)) may overestimate future GPP. 30 We evaluated whether the differences between the observed SCA N BP trends (significant only in L >40N ) and those simulated by the LSMs could be associated with the modelled sensitivities to atmospheric CO 2 concentration (CO 2 ) and growingseason temperature (T) in L >40N (Figure 4). In Figure 4, only models with a too small sensitivity of SCA N BP to T produce a realistic trend of SCA N BP . In contrast, the models indicating sensitivities to T and CO 2 more similar to those estimated by the inversions tend to underestimate the trend in SCA.

Evaluating model biases
Why are the LSM sensitivities of SCA to T positively correlated with their long-term SCA N BP trend (Figure 4), even though CO 2 is a stronger driver of the simulated SCA N BP trend (Figure 3)? We found a clear relationship between the model bias in the trend of SCA N BP and the sensitivity to CO 2 fertilisation in S3 (in line with Wenzel et al. (2016)), but 5 we also found a compensatory effect, where models that overestimate the sensitivity of SCA to T tend to underestimate the sensitivity to CO 2 and vice-versa. LSMs tend to overestimate (underestimate) sensitivity of SCA N BP to T (CO 2 ), compared to the observation-based constraints from inversions. LSMs often compensate too strong (or too weak) simulated water-stress or temperature sensitivity by adjusting photosynthesis parameters (that control CO 2 fertilization) during model optimization to match the observed net terrestrial sink. This compensatory effect has previously been reported by Huntzinger et al. (2012) for 10 the mean terrestrial sink; we find that it could also affect the trends in seasonal CO 2 exchange.
We argue that the trend of SCA N BP can differ between models due to: a) differences in their NPP response to T and CO 2 ; b) differences in turnover times of short-lived C pools by which increased NPP is coupled to increased winter decomposition; (c) phase shifts between GPP and ecosystem respiration. The latter 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 turn-15 over times and associated effective soil depths Koven et al. (2009) (ii) neglect of seasonal acclimation effects on autotrophic and heterotrophic respiration. The sensitivities of NPP to CO 2 and T between models are strongly and consistently correlated with the compensatory effect of the model parameterisations ( Figure S9), but we find no clear relationship between the biases of the modelled SCA N BP trend and the sensitivity of NPP to T (Figure S10), suggesting a key role of respiration. Indeed, the models with SCA N BP trends closer to observations tend to be associated with a lower sensitivity of ecosystem respiration 20 to growing-season temperature ( Figure S4c). Too large turnover of short-lived pools in a model should produce too small increase of the SCA N BP 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 2 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 N BP requires a deeper analysis of transfers between litter and soil 25 organic carbon pools and can be verifiable in future simulations.

Conclusions
Based on our assessment of atmospheric observations and most advanced land-surface model simulations the most likely explanation of the seasonal cycle of atmospheric CO 2 at high latitudes is the CO 2 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 30 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 N BP trend contains valuable information for the turnover times of short-term pools, which await further investigation. Jeong, S.-J., Bloom, A. A., Schimel, D., Sweeney, C., Parazoo, N. C., Medvigy, D., Schaepman-Strub, G., Zheng, C., Schwalm, C. R., Huntzinger, D. N., et al.: Accelerating rates of Arctic carbon cycling revealed by long-term atmospheric CO2 measurements, Science