Dr Xuhui Wang

Dr Vanessa Haverd

Dr Atul Jain

Dr Etsushi Kato

Dr Danica Lombardozzi

Dr Benjamin Poulter

Dr Shilong Piao

Kai Wang

Yilong Wang

Xuhui Wang
email: xuhui.wang@pku.edu.cn

Yue He

Xiangyi Li

Ralph F Keeling

Philippe Ciais

Martin Heimann

Shushi Peng

Frédéric Chevallier

Pierre Friedlingstein

Stephen Sitch

Wolfgang Buermann

Vivek K Arora

Vanessa Haverd

Atul K Jain

Kato

Sebastian Lienert

Danica Lombardozzi

Julia E M S Nabel

Benjamin Poulter

Nicolas Vuichard

Andy Wiltshire

Ning Zeng

Dan Zhu

Shilong Piao

Causes of slowing-down seasonal CO 2 amplitude at Mauna Loa

Keywords: seasonal CO 2 amplitude, slowing down, Mauna Loa, climate change, atmospheric circulation, land

Changing amplitude of the seasonal cycle of atmospheric CO 2 (SCA) in the northern hemisphere is an emerging carbon cycle property. Mauna Loa (MLO) station (20 o N, 156 o W), which has the longest continuous northern hemisphere CO 2 record, shows an increasing SCA before the 1980s (P < 0.01), followed by no significant change thereafter. We analyzed the potential driving factors of SCA slowing-down, with an ensemble of dynamic global vegetation models (DGVMs) coupled with an atmospheric transport model. We found that slowing-down of SCA at MLO is primarily explained by response of net biome productivity (NBP) to climate change, and by changes in atmospheric circulations. Through NBP, climate change increases SCA at MLO before the 1980s and decreases it afterwards. The effect of climate change on the slowing-down of SCA at MLO is mainly exerted by intensified drought stress acting to offset the acceleration driven by CO 2 fertilization. This challenges the view that CO 2 fertilization is the dominant cause of emergent SCA trends at northern sites south of 40°N. The contribution of agricultural intensification on the deceleration of SCA at MLO was elusive according to land-atmosphere CO 2 flux estimated by DGVMs and atmospheric inversions. Our results also show the necessity to adequately account for changing circulation patterns in understanding carbon cycle dynamics observed from atmospheric observations and in using these observations to benchmark DGVMs.

| INTRODUCTION

The seasonal amplitude of atmospheric carbon dioxide (CO 2 ) mole fraction (hereafter SCA) refers to the peak-to-trough magnitude of the detrended seasonal cycle of CO 2 for each calendar year. It results primarily from the convolution of seasonal variations in atmospheric transport and of seasonal variations in terrestrial CO 2 fluxes fixed by photosynthesis and released by respiration and disturbances [START_REF] Bacastow | Seasonal amplitude increase in atmospheric CO₂ concentration at Mauna Loa, Hawaii, 1959-1982[END_REF][START_REF] Keeling | Increased activity of northern vegetation inferred from atmospheric CO₂ measurements[END_REF]. Although the seasonality of ocean CO 2 fluxes and fossil fuel emissions also affect the SCA, their contributions are expected to be much smaller [START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF][START_REF] Keeling | Increased activity of northern vegetation inferred from atmospheric CO₂ measurements[END_REF]. As a result, the SCA is viewed as an integrated signal for studying the response of the terrestrial carbon cycle to global change [START_REF] Buermann | The changing carbon cycle at Mauna Loa Observatory[END_REF][START_REF] Dargaville | Evaluation of terrestrial carbon cycle models with atmospheric CO₂ measurements: Results from transient simulations considering increasing CO₂, climate, and land-use effects[END_REF][START_REF] Forkel | Enhanced seasonal CO 2 exchange caused by amplified plant productivity in northern ecosystems[END_REF][START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF][START_REF] Heimann | Evaluation of terrestrial carbon cycle models through simulations of the seasonal cycle of atmospheric CO₂: First results of a model intercomparison study[END_REF][START_REF] Keeling | Increased activity of northern vegetation inferred from atmospheric CO₂ measurements[END_REF][START_REF] Piao | On the causes of trends in the seasonal amplitude of atmospheric CO₂[END_REF]Zeng et al., 2014).

A global network of more than 20 atmospheric CO 2 monitoring stations was established in the 1980s. During the past three decades, data from this network have shown the SCA increases in the high northern latitudes, but does not increase at mid and low latitude stations [START_REF] Forkel | Enhanced seasonal CO 2 exchange caused by amplified plant productivity in northern ecosystems[END_REF][START_REF] Piao | On the causes of trends in the seasonal amplitude of atmospheric CO₂[END_REF]. A longer perspective on the SCA trend can be gained from observations made at the Mauna Loa station (MLO), located at 19.5°N, 155.6°W and an elevation of 3397 meters in the Pacific Ocean, which has been recording data since 1958 [START_REF] Keeling | Increased activity of northern vegetation inferred from atmospheric CO₂ measurements[END_REF]. Data from MLO show a significant increase in SCA before the 1980s, but this increment then slows down (Figure 1; Figure S1), with the trend after the 1980s being similar to other northern mid-latitude stations [START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF]. Less attention has been paid to this deceleration of SCA at MLO [START_REF] Buermann | The changing carbon cycle at Mauna Loa Observatory[END_REF] than the continuous SCA rise in the Arctic [START_REF] Bacastow | Seasonal amplitude increase in atmospheric CO₂ concentration at Mauna Loa, Hawaii, 1959-1982[END_REF][START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF]. Understanding the deceleration of SCA at MLO is important because the SCA trend is one of the few emergent properties used to constrain future projections of the global carbon cycle (Wenzel, Cox, Eyring, & Friedlingstein, 2016), and the deceleration of SCA at MLO may also offer some indications on whether the persistent increase of SCA at higher latitudes will continue in the future [START_REF] Bacastow | Seasonal amplitude increase in atmospheric CO₂ concentration at Mauna Loa, Hawaii, 1959-1982[END_REF][START_REF] Keeling | Increased activity of northern vegetation inferred from atmospheric CO₂ measurements[END_REF].

In this paper, we simulated the SCA trend at MLO in different time windows during 1959-2016 using an ensemble of land carbon cycle models, ocean carbon cycle models and fossil fuel CO 2 emissions dataset with an atmospheric transport model. We separated the contribution of terrestrial and ocean carbon fluxes, land use emissions, fossil fuel emissions and atmospheric circulation to

Accepted Article

This article is protected by copyright. All rights reserved the SCA trend and trend change. The effect of climate change and CO 2 fertilization on the contribution of terrestrial carbon fluxes were also quantified.

| MATERIALS AND METHODS

| Datasets

| CO 2 data

The weekly in situ carbon dioxide (CO 2 ) mole fraction data at Mauna Loa (MLO) for the period March 1958 to December 2018 were obtained from the Scripps CO 2 Program (http://scrippsco2.ucsd.edu/data/atmospheric_co2/mlo) [START_REF] Keeling | Exchanges of Atmospheric CO₂ and 13 CO₂ with the Terrestrial Biosphere and Oceans from 1978 to 2000. I. Global Aspects[END_REF]. Since calculation of the CO 2 amplitude should only use full year records, only data for the period January 1959 to December 2018 were used in this study. We also used monthly CO 2 dry air mole fractions data for other sites from the NOAA ESRL Carbon Cycle Cooperative Global Air Sampling Network (ftp://aftp.cmdl.noaa.gov/data/trace_gases/co2/flask/surface/) [START_REF] Dlugokencky | Atmospheric Carbon Dioxide Dry Air Mole Fractions from the NOAA ESRL Carbon Cycle Cooperative Global Air Sampling Network[END_REF]. 20 sites with data for at least 25 years in northern latitudes were used.

| Palmer Drought Severity Index (PDSI)

PDSI is a widely used standardized index to describe drought severity. Lower values of PDSI indicate drier conditions [START_REF] Dai | Estimates of Freshwater Discharge from Continents: Latitudinal and Seasonal Variations[END_REF]. We used the Self-calibrating Palmer Drought Severity Index provided by the Climatic Research Unit (CRU), University of East Anglia (https://crudata.uea.ac.uk/cru/data/drought/) [START_REF] Osborn | Monitoring global drought using the self-calibrating Palmer Drought Severity Index [in" State of the Climate in 2016[END_REF][START_REF] Van Der Schrier | A scPDSI-based global data set of dry and wet spells for 1901-2009[END_REF].

Monthly PDSI for the period 1901 to 2016 was calculated from the CRU TS 3.26 monthly climate data which has a spatial resolution of 0.5º×0.5º.

| Dynamic global vegetation models (DGVMs)

Land-atmosphere CO 2 exchange (net biome productivity, NBP) was provided by 12 dynamic global vegetation models (DGVMs) (Table S1). All the model simulations were performed following the TRENDY inter-comparison protocol for the common period of 1901-2016 (TRENDYv6) using the same climate drivers obtained from CRU-NCEP v8 [START_REF] Harris | Updated high-resolution grids of monthly climatic observations -the CRU TS3.10 Dataset[END_REF]Wei et al., 2014), atmospheric CO 2 values from a combination of ice core records and atmospheric observations, and land use change from the HYDE data set (Le Quéré et al., 2018). Simulations were performed for three scenarios (S1, S2, S3) to facilitate isolating the

Accepted Article

This article is protected by copyright. All rights reserved effects of increasing CO 2 mole fraction, climate change, and land use [START_REF] Sitch | Recent trends and drivers of regional sources and sinks of carbon dioxide[END_REF]. In S1, only atmospheric CO 2 mole fraction changed with time, while climatology without any trend was used as climate forcing and other factors were kept constant. In S2, both atmospheric CO 2 mole fraction and climate varied over time, while S3 considered varying atmospheric CO 2 mole fraction, climate, and land use. The varying climate accounts for the historical change of radiation, precipitation, surface air temperature, surface pressure, air humidity and wind speed. The 12 DGVMs use some (depending on the requirement of model inputs) of these climate variables to drive the simulation of NBP. Consequently, the effects of CO 2 , climate, and land use could be estimated from S1, S2-S1 and S3-S2, respectively. Monthly net biome productivity (NBP), gross primary productivity (GPP) and total ecosystem respiration (TER) in S1, S2, S3 and soil moisture in S2 were used in this study. Note that we used the soil moisture output from all of the models except for CLASS-CTEM due to the missing soil moisture output. Climate forcing of the DGVMs was provided by CRU-NCEP v8, which is a globally gridded (0.5º×0.5º) data set based on meteorological observations and climate reanalyses [START_REF] Harris | Updated high-resolution grids of monthly climatic observations -the CRU TS3.10 Dataset[END_REF]Wei et al., 2014), for the period 1901-2016. This data set also provided values of monthly mean temperature, total rainfall and downward shortwave radiation for the period May to September for the analyses of the correlation between interannual variations of climate variables and carbon fluxes (NBP, GPP and TER).

| Fossil fuel CO 2 emissions

The gridded monthly fossil fuel CO 2 emission data from the Carbon Dioxide Information Analysis Center (CDIAC) spanning the period between 1959 and 2013 was used in this study, as it is also used by the Global Carbon Project for estimating fossil fuel emissions. The fossil fuel emissions include emissions from solid, liquid and gas fuels, and the seasonal cycles were estimated by using a proportional-proxy method [START_REF] Andres | Monthly, global emissions of carbon dioxide from fossil fuel consumption[END_REF]. For 21 countries (covering about 80% of total global emissions), fossil fuel CO 2 emission data were obtained from available monthly fossil fuel consumption data. For the other countries, for which such monthly data were not available, the seasonal cycle of fossil fuel emissions were represented by the available data in view of climatic and economical similarities (also geographic proximity for several countries). For some years there were gaps in the data, and to fill these, Monte Carlo methods were used to create values from years with known monthly fractions [START_REF] Andres | Monthly, global emissions of carbon dioxide from fossil fuel consumption[END_REF]. Since the CDIAC data were only available until 2013, they were extrapolated to 2016 by

Accepted Article

This article is protected by copyright. All rights reserved using global fuel consumption data provided by the Open-source Data Inventory for Anthropogenic CO 2 (ODIAC) for 2014-2016 [START_REF] Oda | The Open-source Data Inventory for Anthropogenic CO₂, version 2016 (ODIAC2016): A global monthly fossil fuel CO₂gridded emissions data product for tracer transport simulations and surface flux inversions[END_REF]. To assess the impacts of uncertainties associated with estimates of fossil fuel emissions on the SCA trend change, we also used two alternative gridded fossil fuel CO 2 emission datasets: CEDS [START_REF] Hoesly | Historical (1750-2014) anthropogenic emissions of reactive gases and aerosols from the Community Emissions Data System (CEDS)[END_REF] and PKU [START_REF] Wang | High-resolution mapping of combustion processes and implications for CO 2 emissions[END_REF].

| Ocean-atmosphere CO 2 exchange

The ocean-atmospheric CO 2 exchange was simulated by the PlankTOM5 model combined with the NEMO model [START_REF] Buitenhuis | Biogeochemical fluxes through microzooplankton[END_REF]. PlankTOM5 is a biogeochemical model forced by compounds in rivers, sediments and dust [START_REF] Aumont | An ecosystem model of the global ocean including Fe, Si, P colimitations[END_REF][START_REF] Da Cunha | Potential impact of changes in river nutrient supply on global ocean biogeochemistry[END_REF]; NEMO is a global ocean general circulation model, forced by daily wind and precipitation data from the National Center for Environmental Prediction (NCEP) reanalysis [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF]. To assess the impacts of uncertainties associated with ocean fluxes on the SCA trend change, we also used six other ocean models from the global carbon project (Le Quéré et al., 2018), namely: CCSM-BEC, CSIRO, MPIOM, NEMO-PISCES_CNRM, NorESM and RECOM.

| Global atmospheric transport model

We estimated the atmospheric CO 2 mole fraction at MLO by forcing the global atmospheric transport model of the Laboratoire de Météorologie Dynamique version 5 [START_REF] Hourdin | The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection[END_REF]Locatelli et al., 2015) (LMDZv5) with the inputs of land-atmosphere CO 2 exchange, fossil fuel CO 2 emissions and ocean-atmosphere CO 2 exchange fluxes. The LMDZ transport model was run on a 1.875º×3.75º (latitude × longitude) horizontal grid, with 39 layers between the surface and the top of the atmosphere. LMDZ was nudged to European Centre for Medium-Range Weather Forecasts (ECMWF) reanalyzed winds, and was run in an off-line mode in which transport mass fluxes were read from pre-computed archives from the same model rather than computed online.

The physical convection scheme used in this study is also used in the LMDZ5A model in the CMIP5 experiments. Note that the time span of each ECMWF wind product did not cover the whole study period of 1959-2016. Therefore we used two ECMWF reanalysis products covering different periods, i.e. ERA-20C for 1959-2010, and ERA-Interim for 1979-2016 (https://www.ecmwf.int/en/forecasts/datasets). The ERA-20C and ERA-Interim products differ in

Accepted Article

This article is protected by copyright. All rights reserved their spatial resolution, model orography, assimilated data, assimilation methodology, etc., so that the wind fields in ERA-20C and ERA-Interim are also different [START_REF] Poli | ERA-20C: An Atmospheric Reanalysis of the Twentieth Century[END_REF]. We compared the modeled CO 2 seasonal amplitudes during 1980-2010 from simulations using ERA-20C (MEC) and ERA-Interim (MEI) wind data sets. Figure S2 shows that the modeled CO 2 seasonal amplitudes are very close to each other and are not sensitive to the choice of the ECMWF product.

Therefore, for each transport simulation (see below for details), we combined the modeled CO 2 seasonal amplitudes using ERA-20C for 1959-1979and those using ERA-Interim for 1980-2016 to form a full time series of modeled SCA from 1959 to 2016. The analysis of the southern boundary of the Hadley cell shows no drastic change in the combination of the two-reanalysis products (Figure S3).

We were careful to ensure that the use of the two different wind field products did not lead to a trend change of SCA. Firstly, we verified that the SCA values calculated with the different wind fields were consistent for the overlapping time period when both wind field data sets were available (Figure S2). Secondly, we found that the change in the trend of SCA was also detected when the ERA-20C transport was used (Figure S4). Moreover, when we changed the start year for using ERA-Interim to a year in the late 1980s or early 1990s, such as 1990, we found that the SCA trend change still occurred in the 1980s (Figure S4). The matrices of SCA trend change (the orange patches in Figure S4) did not move to 1990 to match the change in the ERA-Interim start year. These results indicate that the contribution of the change in atmospheric circulation to the SCA trend change is not induced by use of the two wind field products to generate the single SCA time series.

To estimate the robustness of the trend change of SCA induced by the atmospheric circulation, we used reanalyzed wind fields, from both ECMWF and NCEP [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF][START_REF] Kistler | The NCEP-NCAR 50-Year Reanalysis: Monthly Means CD-ROM and Documentation[END_REF] in the transport simulations. However, the Simulation T1 used time-varying fluxes and NCEP R1 wind field forcing (T1 NCEP ), simulated a SCA trend of 0.01 ppm yr -1 (P = 0.05), largely underestimating the observed SCA increase of 0.04 ppm yr -1 (P < 0.01) during 1959-1984 (Figure S5a), and much lower than the simulated SCA trend from the Simulation T1 (0.03 ppm yr -1 , P < 0.01) forced with the ECMWF wind field (T1). This result illustrates the fact that circulation changes in the data-sparse period before the 1980s are more uncertain than in recent decades [START_REF] Fujiwara | Introduction to the SPARC Reanalysis Intercomparison Project (S-RIP) and overview of the reanalysis systems[END_REF] and the SCA trend seems to be highly sensitive to the atmospheric

Accepted Article

This article is protected by copyright. All rights reserved circulation changes. The differences in the wind field between ECMWF and NCEP R1 before the 1980s may come from many sources, such as forecast models, input data sets and data assimilation systems for example. Although it is difficult and beyond our scope to evaluate the reanalyzed wind fields, we hypothesize that differences in the input sea surface temperature (SST) could be an important driver for wind field differences between the two reanalyses: The HadISST used by ECMWF represents a major improvement over its predecessor GISSST, which is used by NCEP R1, in forcing tropical atmospheric circulations [START_REF] Rayner | Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century[END_REF]. Differences between HadISST and NCEP R1 are mainly found before the 1980s [START_REF] Rayner | Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century[END_REF]. Since the observed change in the SCA trend can only be reproduced by the simulation using the wind products from ECMWF (T1), and not by the simulation using NCEP R1 (T1 NCEP ), the analysis of the impact of atmospheric circulation in this study is based only on the ECMWF data. We also simulated CO 2 mole fraction to estimate SCA for other 20 sites coving a range of northern latitudes.

| Analysis approach

| Design of transport experiments

To separate the effects of six factors that can affect the CO 2 seasonal amplitude trends: changes of atmospheric CO 2 ('CO 2 '), climate ('CLIM'), land use ('LU'), fossil fuel ('FF'), ocean carbon flux ('Ocean') and atmospheric circulation ('Wind'), seven transport simulations (T1-T7) were designed, using global transport models coupled with DGVMs (Table 1). The first (T1) used time-varying S3 monthly NBP, fossil fuel CO 2 emissions, and ocean-atmosphere CO 2 exchange, coupled with the use of time-varying wind fields in the LMDZv5. The T1 simulation thus integrates the effects of all six factors. The second simulation (T2) is driven by time-varying S1 monthly NBP and constant fossil fuel CO 2 emissions and ocean-atmosphere CO 2 exchange, coupled with the LMDZv5 transport model with variable winds, which reflects the combined effects of 'CO 2 ' and 'Wind'. The third simulation (T3) is driven by time-varying S2 monthly NBP, while the other factors remain the same as in T2. T4 is driven by time-varying S3 monthly NBP, with the other factors kept the same as in T3. Simulation T5 is forced by historical varying wind, but uses constant CO 2 fluxes (land-atmosphere CO 2 exchange, fossil fuel CO 2 emissions, and ocean-atmosphere CO 2 exchange) with values from 1959. In Simulations T6 and T7, only fossil fuel CO 2 emissions and ocean-atmosphere CO 2 exchange are varying respectively, both coupled with the varying winds in LMDZv5. Therefore, the difference between two specific

Accepted Article

simulations indicates the effect of one factor. That is, the effects of these factors atmospheric CO 2 ('CO 2 '), climate ('CLIM'), land use ('LU'), fossil fuel ('FF') and ocean carbon flux ('Ocean'), can be calculated by T2-T5, T3-T2, T4-T3, T6-T5 and T7-T5, respectively. The effect of atmospheric transport ('Wind') is simply obtained from T5.

To investigate the impact of fossil fuel emissions on the change in the trend of SCA, we run additional T5 and T6 simulations with three gridded fossil fuel CO 2 emission maps: CDIAC (used in the reference simulation), CEDS and PKU. In these additional simulations, we only used LPX for the land NBP fluxes. For the ocean, we used 7 models from the global carbon project, namely: CCSM-BEC, CSIRO, MPIOM, NEMO-PISCES_CNRM, NEMO-PlankTOM5 (used in the reference simulation), NorESM and RECOM. We ran additional T5 and T7 simulations using different ocean fluxes.

To investigate the impact of agriculture on the SCA at MLO, we conducted additional experiments by using a cropland mask for the year 2014 from MODIS (MCD12C1 v006, https://lpdaac.usgs.gov/products/mcd12c1v006/). In each grid box of the transport model, if the fraction of cropland is larger than 30% (in most grid boxes, cropland is the dominant land cover type), the grid box is treated as "cropland". For the 12 DGVMs, we run additional simulations (TC) where the fluxes in the cropland grid boxes cycle with the 1959 monthly fluxes, while the fluxes in the other grid boxes are the same as for Simulation T1. Therefore, the modeled SCA trend from TC simulation is not contributed by the cropland, that is, agricultural flux. The difference between the T1 and TC simulations is, the DGVMs' simulated impact of agriculture. In addition, we also investigate the impact of agriculture, based on inversion results, namely: Jena_s57X (Rödenbeck, Zaehle, Keeling, & Heimann, 2018a), Jena sEXT (Rödenbeck, Zaehle, Keeling, & Heimann, 2018a), Jena s85 [START_REF] Rödenbeck | How does the terrestrial carbon exchange respond to inter-Annual climatic variations? A quantification based on atmospheric CO 2 data[END_REF] and CAMS (https://apps.ecmwf.int/datasets/data/cams-ghg-inversions/) using the same cropland mask.

Jena_s57X and Jena_sEXT cover the period from 1957 to the present day, while Jena_s85 and CAMS only start in 1985 and 1979 respectively. We run the transport model, both with the actual fluxes provided by the inversions (TI), and with the fluxes in the cropland grid boxes cycling through the 1959 monthly fluxes while the fluxes in the other grid boxes are the same as in Simulation TI (TIC). Since atmospheric inversions were driven by atmospheric CO 2 observations, which convoluted all land-atmosphere CO 2 exchange, if agricultural intensification contributes to

Accepted Article

This article is protected by copyright. All rights reserved SCA trend, it should be manifested by the contribution of NBP over agricultural zone to SCA change, represented by difference between TI and TIC.

| Estimation of seasonal CO 2 amplitude (SCA) and SCA trend

SCA derived from the MLO observatory and simulated data was calculated using a standard tool, CCGCRV [START_REF] Thoning | Atmospheric carbon dioxide at Mauna Loa Observatory: 2. Analysis of the NOAA GMCC data, 1974-1985[END_REF]). The long term curves and annual oscillation of atmospheric CO 2 were firstly obtained by fitting a function to the raw data. The function consists of a quadratic polynomial for the long-term trend, and four-harmonics for the annual cycle. Then a 50-day cutoff value short term filter and a 667-day cutoff value long term filter were applied to the residuals between the raw data and the fitted function, so that any short term variations or annual cycles that are still present in the residuals after fitting the function are removed [START_REF] Thoning | Atmospheric carbon dioxide at Mauna Loa Observatory: 2. Analysis of the NOAA GMCC data, 1974-1985[END_REF]. The detrended seasonal cycle is obtained by adding the filtered residuals using the short term cutoff value to the annual cycle. The peak-to-trough amplitude was calculated as the difference between the maximum and the minimum CO 2 mole fraction of the CO 2 seasonal cycle in each calendar year. In general, the maximum and minimum of the CO 2 seasonal cycle at MLO appeared in May and at the end of September or early October, respectively, which is consistent with the results of previous studies [START_REF] Buermann | The changing carbon cycle at Mauna Loa Observatory[END_REF][START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF]. As we focus on the change of SCA, the anomalies of SCA rather than actual values were used in the analysis.

To ensure that the observed slow-down in the increasing trend of SCA was not an artifact introduced in some way by the specific tool used to extract the seasonal trend, we repeated the SCA time series extraction process using two other algorithms: Seasonal Trend decomposition using LOESS (STL, where LOESS is an abbreviation for locally weighted scatterplot smoothing, available at https://stat.ethz.ch/R-manual/R-devel/library/stats/html/stl.html) and HPspline [START_REF] Pickers | Investigating bias in the application of curve fitting programs to atmospheric time series[END_REF]. STL uses a moving average algorithm to fit the seasonal cycle, whereas the HPspline algorithm involves fitting data to a harmonic function, a polynomial equation, and a stiff cubic spline. The full details of the three algorithms, and of their different performances in extracting long-term trends and seasonal cycles, are discussed in the previous study [START_REF] Pickers | Investigating bias in the application of curve fitting programs to atmospheric time series[END_REF]. The use of STL or HPspline to extract the SCA time series did not lead to any change in the observed deceleration of the SCA at MLO (Figure S6).

We tested the impacts of different breakpoints and anomalous years on the estimation of SCA

Accepted Article

This article is protected by copyright. All rights reserved trend. Different breakpoints from 1983 to 1986 and exclusion of anomalous years were considered. Anomalous years are identified as El Niño years (i.e., 1983-1984, 1987-1988, 1997[START_REF] Heimann | Evaluation of terrestrial carbon cycle models through simulations of the seasonal cycle of atmospheric CO₂: First results of a model intercomparison study[END_REF]and 2015-2016), La Niña years (i.e., 1970-1971, 1974-1975and 2010-2011) and years when there are large volcanic eruptions (i.e., 1963, 1980, 1982, 1991 and 2011). In this sensitivity test, we remove SCA of all possible combinations of anomalous years (from any single year to all the years) in trend derivation (Figure 2). The frequency distribution of trends is obtained from 5000 bootstrap analyses. In addition, we fit linear trends to the SCA at MLO, using at least a 20-year time span, with different combinations of start years (ranging from 1959 through 1997), and end years (ranging from 1978 through 2016) (Figure 3).

| CO 2 source region (footprint) analysis

To identify the CO 2 flux regions that affect the change of SCA at MLO, we calculated the sensitivity of the monthly mean CO 2 mole fraction in September at MLO to land-atmosphere/ocean-atmosphere CO 2 exchange between May and September (footprint) by using the adjoint code of LMDZv5 [START_REF] Chevallier | Inferring CO₂ sources and sinks from satellite observations: Method and application to TOVS data[END_REF]. Because the seasonal amplitude is computed as the difference between the monthly mean CO 2 mole fraction in September and that in May, which is a result of the cumulative sink of the terrestrial and oceanic biosphere during this period, this footprint represents the amount by which the monthly mean CO 2 mole fraction will drop from the May CO 2 mole fraction, when the land-atmospheric/ocean-atmospheric CO 2 exchange decreases by 1 kg C m -2 h -1 (here decrease means a larger sink in the biosphere) for each day from May to September. The adjoint code analytically computes the partial derivatives of all processes within LMDZv5 following the chain rule and thus allows us to obtain partial derivatives of monthly mean CO 2 mole fraction with respect to the land-atmosphere/ocean-atmosphere CO 2 exchange for all grid points. The footprints in units of mole fraction (ppm) per unit flux (kg C m -2 h -1 ) were computed every year for all global grid-point fluxes of the transport model at the daily scale since May 1 st . We then summed the footprints for all the days from May 1 st to September 30 th in each year between 1959 and 2016 to get a total footprint of the flux regions affecting the CO 2 drawdown, and hence the SCA, at MLO from May to September. The main CO 2 source regions estimated using LMDZv5 are similar to those found by using TM2 [START_REF] Kaminski | Sensitivity of the seasonal cycle of CO₂ at remote monitoring stations with respect to seasonal surface exchange fluxes determined with the adjoint of an atmospheric transport model[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved To examine whether the identified terrestrial CO 2 source regions are indeed the regions mostly responsible for the trend change of SCA at MLO, we perform a separate transport experiment: Simulation T8. In T8, we keep the NBP of lightly contributing regions (those with mean surface flux sensitivity less than 150 ppm per kg C m -2 h -1 ) constant at their 1959 values throughout the 1960-2016 simulation period, while allowing the NBP of the more heavily contributing regions (with mean surface flux sensitivity larger than 150 ppm per kg C m -2 h -1 ) to vary according to the NBP simulated by the 12 DGVMs. By comparing T8 and T2, we found that the identified source regions almost fully explain the SCA trend change at MLO between 1959MLO between -1984MLO between and 1985MLO between -2016. . On the other hand, in some regions, the reversal of trend of NBP driven by climate change To evaluate whether the oceanic contribution to seasonal CO 2 mole fraction at MLO has changed, due to changes in the atmospheric circulation during the period 1959-2016, we calculated the area weighted sum of sensitivity of CO 2 at MLO to the ocean flux (ocean footprint) and that to the land flux (land footprint), respectively. The ratio of the ocean footprint to the land footprint is calculated for the period 1959-2016 for the following regions: the entire globe; the footprint regions where the mean surface flux sensitivity was larger than 150 ppm per kg C m -2 h -1 ; and for each latitudinal band.

Accepted Article

Model weights in the averaging

The model average was derived using two approaches: arithmetic model averaging (AMA) and Bayesian model averaging (BMA). In AMA, each model was given equal weight, without considering the ability of each model to accurately reproduce the observations. In BMA, in contrast, the weight assigned to each model (listed in table S1) was determined by the ability of the model to reproduce the observations, with the sum of the weights being equal to one. In our case, models (Simulation T1), having higher correlations with the observations, were assigned with higher weights in the BMA scheme. The better the agreement of the amplitude anomalies of a model with the observed data, the greater the weight assigned. The optimal probability density function was calculated by Monte Carlo Markov chains [START_REF] Vrugt | Ensemble Bayesian model averaging using Markov Chain Monte Carlo sampling[END_REF]. Three models (LPX, ISAM and CABLE) had high weights up to 90% (table S1). To investigate the robustness of BMA, weights of simulated SCA using BMA were derived from other independent observations: the interannual variability of CO 2 amplitude and the mean CO 2 seasonal cycle. Simulated SCA trend using BMA and the correlation between the observed SCA and simulated SCA were estimated. Uncertainties of the BMA-based SCA trends were estimated by 5000 bootstrap estimates. First, a bootstrap sample set of size 1000 was composed of SCA trends from 12 individual DGVMs. The number of samples for each model is based on its BMA weight. For example, the weight of CABLE is 0.12 (table S1), so 120 samples were taken with the SCA trend estimated by the CABLE model. Then, we resample this bootstrap sample set 5000 times and calculate the standard deviation of those bootstrap estimates. Note that when the trend of one model during 1959-1984 was sampled, the trend of the same model during 1985-2016 was sampled simultaneously.

| RESULTS

| Observed and simulated SCA trend at MLO

SCA at MLO increased by 1.2 ppm over past six decades (Figure 1), but the increasing trend is not homogeneous. The increasing trend was 0.04±0.01 ppm yr -1 (P < 0.01) from 1959 to the mid-1980s, but reduced by 50% to 0.02±0.01 ppm yr -1 (P > 0.05) from the mid-1980s to 2018 (Figure 1). To assess the robustness of the slow-down in SCA trend at MLO since the 1980s, we test the histograms of SCA trend before and after the 1980s considering different breakpoints (Figure 2; Figure S1). Although trend estimates might be affected by few extreme values, for

Accepted Article

This article is protected by copyright. All rights reserved example, the historically high values of SCA in 2015 and 2016 (Figure 1) [START_REF] Bastos | Impact of the 2015/2016 El Niño on the terrestrial carbon cycle constrained by bottom-up and top-down approaches[END_REF], whether including these two extreme years does not change the conclusion that SCA trend at MLO slows down (Figure 2; Figure S1). This also holds true when SCA of the recent two years (2017 and 2018), whose SCA drops down below that of 2014, were excluded (Figure 1; Figure 2). In order to systematically estimate the impact of extreme values on trend estimates, we perform bootstrap analyses by randomly excluding anomalous years in deriving the trend. The SCA trend significantly slows down by 0.02±0.01 ppm yr -1 (P < 0.01) after the mid-1980s (Figure 2). We also fit linear curves to the SCA at MLO, using windows with at least 20 years and as shown in Figure 3a, the SCA trend during the first three decades is significantly positive (> 0.03 ppm yr -1 , P < 0.01), but, after the 1980s, the trends become insignificantly different from zero, or even negative. In addition, three different methods to extract SCA from observed CO 2 mole fraction at MLO (frequency based time series decomposition (CCGCRV), seasonal trend decomposition using locally weighted scatterplot smoothing (STL), and time series decomposition combining harmonic functions, polynomial and cubic spline (HPspline)) show similar slow-down of SCA trend at MLO since the mid-1980s (Figure S6). We conclude that the slow-down in the SCA trend at MLO is robust to potential artifacts induced by breakpoints, anomalous events and methods to extract SCA.

To understand the mechanisms behind the deceleration of the SCA at MLO, we compiled gridded estimates of ocean CO 2 fluxes, terrestrial CO 2 flux (or net biome productivity, NBP) and fossil fuel emissions, and computed the pertaining CO 2 mole fraction using the LMDZv5 tracer transport model [START_REF] Chevallier | Inferring CO₂ sources and sinks from satellite observations: Method and application to TOVS data[END_REF]. Uncertainties in the spatiotemporal patterns of NBP, which control SCA, are addressed by using monthly NBP from the TRENDYv6 dynamic global vegetation models (DGVMs) (Table S1) (Le Quéré et al., 2018). This 12-model NBP ensemble contains factorial experiments executed by each DGVM. Seven transport model simulations (T1-T7) were performed to isolate the contributions to SCA trends of the driving factors: terrestrial NBP (changes in atmospheric CO 2 , climate and land use), fossil fuel emissions, ocean CO 2 fluxes, and atmospheric circulation (Table 1, Table S2).

First, we examined whether SCA simulated using NBP derived from all factors (Simulation T1) can reproduce the observed SCA trend. To attribute the observed trend change, we used model

Accepted Article

This article is protected by copyright. All rights reserved ensemble output calculated using Bayesian Model Averaging (BMA), which has been proved to result in more robust prediction than the arithmetic mean of the model ensemble (AMA). Unlike simple model averaging, our BMA ensemble gives more weight to those models better reproducing the observed variations in SCA (Figure 1 and Figure S7). The BMA model ensemble of the 12 models can reproduce the deceleration of SCA since the 1980s (Figure 3b), and is also able to reproduce the latitudinal gradients of the observed SCA trend (Figure S8). The model weights in the BMA ensemble can also be derived based on the model's performance in reproducing either the detrended interannual variations of SCA or the mean seasonal cycle of atmospheric CO 2 at MLO. In these two alternative methods of weight-deriving, using different characters of SCA, the alternative BMA ensembles perform similarly to the original BMA ensemble (Figure S9), illustrating the robustness of using BMA to hindcast and attribute the SCA change. Therefore, the BMA of the model ensemble was used to explore the contribution of each potential driving factor to the deceleration of SCA based on factorial transport simulations (Table 1, Table S2).

| Drivers of SCA slowing-down

| Climate change

The impact of the response of NBP to climate change on SCA was derived from the difference between Simulation T3 and Simulation T2 (Table 1, Table S2). The BMA of the DGVMs shows that the trend in NBP response to climate change drove SCA to increase significantly (> 0.01 ppm yr -1 , P < 0.05) from the 1960s to the 1980s, and that its effect then weakened or reversed (~ 0 ppm yr -1 , P > 0.10) (Figure 3c). This decrease of SCA trend driven by climate change is consistent with the observations (Figure 3). When examining the performance of individual DGVMs, we found that the three models that are most successful in reproducing the variations of SCA (R 2 > 0.31), i.e., the ones with larger weights in the BMA ensemble (Table S1), consistently show a stall or decrease of climate-driven SCA since the 1980s (Figure S10).

To diagnose regions, and possible mechanisms, responsible for this climate-induced deceleration of SCA, we examine the change of NBP driven by climate change (NBP clim ), over the land areas which play the most important role in the SCA at MLO: East Asia and western North America (Figure 4a,b). We found a significant reversal of the NBP clim trend over western North America and Eastern China (Figure 4d,e; Figure S11b). By performing another transport experiment

Accepted Article

This article is protected by copyright. All rights reserved (Simulation T9, Table 1), we found that the NBP clim trend reversal over these regions, which is particularly obvious in East Asia and western North America, can almost fully explain the climate change-induced slow-down in the SCA trend at MLO since the 1980s (Figure S11). Furthermore, we found that the reversal of the NBP clim trend over western North America is consistent with the aggravated drought stress on NBP there during the past two decades [START_REF] Méndez | Regional aspects of prolonged meteorological droughts over Mexico and central America[END_REF][START_REF] Schwalm | Reduction in carbon uptake during turn of the century drought in western North America[END_REF]. Indeed, trends in the Palmer Drought Severity Index (PDSI, low values indicate drier conditions) [START_REF] Osborn | Monitoring global drought using the self-calibrating Palmer Drought Severity Index [in" State of the Climate in 2016[END_REF][START_REF] Van Der Schrier | A scPDSI-based global data set of dry and wet spells for 1901-2009[END_REF] and simulated soil moisture changed from increasing to decreasing over the western USA after the 1980s (Figures 5). Over most of Eastern and Northeast China, changes in PDSI and simulated soil moisture also consistently show intensified drought during the last three decades (Figures 5), in parallel with a faster warming trend and a stall or decrease in precipitation since the 1980s (Figure S12). Despite relatively low surface flux sensitivity in South America, the significant reversal of NBP trend in South America may also contribute the deceleration of SCA at MLO (Figure 4). Using a sensitivity simulation (T9 S.Am ), we confirm the sizeable contribution (-0.01±0.01 ppm yr -1 ) of regions where NBP clim trend reversed since the 1980s in South America to the decrease of SCA trend at MLO, which also coincides with the intensified drought stress since the 1980s (Figure 5). As the activity of vegetation during the growing season can be largely influenced by moisture availability (Figure S13), such increasing drought stress could have contributed to the decrease of photosynthetic uptake (gross primary productivity, GPP) since the 1980s (Figure S14), and thus of NBP over the footprint region of MLO (Figure 4f).

| Change in atmospheric circulations

In addition to NBP response to climate change, we find that changes in atmospheric circulation could also have contributed to the trend change of SCA at MLO by shifting the origin of the CO 2 . This mechanism is investigated by using LMDZv5 transport model simulations with an inter-annual varying wind field and constant CO 2 flux from 1959 (Simulation T5, Table 1). From the 1960s to the 1980s, the changing atmospheric circulation contributed to the trends in SCA of between 0.01 ppm yr -1 and 0.02 ppm yr -1 (Figure 3h). The change was due to a northward movement of the regions influencing MLO in the early 1980s (Figure S15). However, after the 1980s, the change in atmospheric circulation led to a significant decrease of the SCA at MLO (~ -0.01 ppm yr -1 ) (Figure 3h; Figure S15a). During 1959During -1984, atmospheric , atmospheric

circulation change

Accepted Article

This article is protected by copyright. All rights reserved contributed by 0.01 ppm yr -1 to the SCA trend, while contributed by -0.01 ppm yr -1 during 1985-2016. Therefore, the impact of change in atmospheric circulation acts as strong as that of climate change on trend change of SCA at MLO. The negative contribution of the decadal variations in atmospheric circulation is robust to patterns of NBP simulated by DGVMs, since 8 of the 12 models, accounting for more than 95% of the weight in the BMA model ensemble, consistently suggest a circulation-induced reduction in SCA (Figure S16). Two processes are potentially responsible for the circulation-induced change in SCA. The first is the changing relative contribution of oceanic and terrestrial carbon fluxes to the CO 2 mole fraction at MLO. We found that the contribution of land carbon fluxes to the CO 2 mole fraction at MLO, particularly that between 30°N and 50°N, increased before the 1980s (Figure S17), which would contribute to the increasing SCA at MLO, since a larger contribution from northern temperate terrestrial carbon fluxes induces a larger seasonal amplitude [START_REF] Lintner | Seasonal circulation and Mauna Loa CO₂ variability[END_REF][START_REF] Piao | Interannual variation of terrestrial carbon cycle: Issues and perspectives[END_REF]. The second process is the north-south movement of the CO 2 source regions of the MLO.

To explore this process, we plotted the year-to-year variations in the wind field (Figure S15a). The results show a noticeable southward movement of the CO 2 source region (footprint) of MLO since the 1980s (Figure 4c; Figure S15a), probably caused by the expansion of the Hadley cell (Figure S3), making it less affected by higher northern latitudes. Since lower latitude lands have less seasonality in NBP (Figure S15b), the southward movement in the source regions, bringing more air mass from lower latitudes, tends to decrease the SCA at MLO, an effect which has been neglected by previous studies. In addition to the wind fields analyzed by ECMWF, we also explored the contribution of circulation changes using the wind fields from the NCEP R1 reanalysis [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF][START_REF] Kistler | The NCEP-NCAR 50-Year Reanalysis: Monthly Means CD-ROM and Documentation[END_REF]. Despite uncertainties in the wind field before the 1980s, the circulation-induced trend of SCA since the mid-1980s is similar with different wind fields (Figure S5), partly supporting the wind field is indeed a driver to deceleration of SCA since the mid-1980s.

| CO 2 fertilization effects

CO 2 fertilization in the increase of SCA can be used to project future photosynthesis (Wenzel, Cox, Eyring, & Friedlingstein, 2016). Increasing CO 2 mole fraction in the atmosphere was found to be a dominant driver for the increase of SCA at northern hemisphere sites south of 40°N [START_REF] Forkel | Enhanced seasonal CO 2 exchange caused by amplified plant productivity in northern ecosystems[END_REF]. We also found a contribution of increasing CO 2 mole fraction to the SCA trend at MLO. The contribution of increasing CO 2 mole fraction did not lead to a deceleration of

Accepted Article

This article is protected by copyright. All rights reserved the SCA, but instead slightly increased the positive trend of SCA over the past three decades (trend increased by 0.01 ppm yr -1 during 1985-2016 compared to that during 1959-1984) (Figure 3d). The accelerated increase of CO 2 mole fraction enhanced the GPP but less affected the seasonality of respiration [START_REF] Keenan | Recent pause in the growth rate of atmospheric CO₂ due to enhanced terrestrial carbon uptake[END_REF], resulting in more cumulative NBP during the net carbon uptake period and the increase of SCA trend since the 1980s.

| Land use change, ocean carbon fluxes and fossil fuel emission

Based on the BMA ensemble of the DGVMs, we found that land use change makes no significant contribution to decadal variations of SCA at MLO (Figure 3e), with its magnitude of contribution is less than 0.001 ppm yr -1 . Previous studies [START_REF] Gray | Direct human influence on atmospheric CO₂ seasonality from increased cropland productivity[END_REF]Zeng et al., 2014) indicate that agricultural intensification, as a component of land use change, was a key driver of the SCA increase at MLO. Due to under-representation of cropland managements in current DGVMs, the DGVM simulated agricultural contribution to SCA trend at MLO is largely uncertain (Zhao et al., 2016). Here, we further performed two additional transport experiments (TI and TIC) with NBP from long-term atmospheric inversion models to quantify the contribution of agriculture zones to SCA trend change at MLO (Table 1, Table S2). The results indicate that the contribution of agriculture zones to SCA trend at MLO has ranged between 0.006 ppm yr -1 and 0.008 ppm yr -1 during 1959-1984 (16-22% of the observed SCA trend), which did not change significantly during 1985-2014 (between 0.004 ppm yr -1 and 0.007 ppm yr -1 , 25-44% of observed SCA trend) (Figure 6). Thus, with the independent evidence from atmospheric inversions, we confirm the role of agricultural intensification in the increase of SCA, but it has little contribution (<0.002 ppm yr -1 ) to the deceleration of SCA at MLO (Figure 6).

The contributions from fossil fuel emissions and ocean carbon fluxes to SCA change at MLO is also not significant (Figure 3f,g; Figures S18,S19). The magnitude of the contributions from fossil fuel emissions to the SCA trend before 1984 is less than 0.001 ppm yr -1 . The contribution from the ocean flux is less than 0.002 ppm yr -1 . After the mid-1980s, land use change and ocean fluxes still contributed little to the SCA trend at MLO. Fossil fuel emissions could lead to a slightly larger increasing trend of SCA (0.004 ppm yr -1 ) after the mid-1980s than before.

| DISCUSSIONS

Accepted Article

This article is protected by copyright. All rights reserved

The increase of SCA at MLO and BRW over the past six decades is a key fingerprint of the current perturbation of land carbon fluxes [START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF], however, the increase of SCA at MLO is not persistent. The differences in trends of SCA at low-latitude MLO (19.5°N) and high-latitude BRW (71.3 °N) over the recent 2-3 decades imply that there are emerging regional differences in ecosystem responses to climate change. Climate change has negative impact on NBP of the temperate northern hemisphere due to intensified drought stress, which induced by an increased vapor pressure deficit associated with the faster warming trend and a stall or decrease in water supply from precipitation. Either the lack of drought stress [START_REF] Nemani | Climate-driven increases in global terrestrial net primary production from 1982 to 1999[END_REF], or the lack of drought stress intensification, in the higher northern latitudes can help to explain why the slow-down of the SCA trend took place in MLO but not in Point Barrow (BRW,71.3 °N), where the trend of seasonal amplitude was mainly contributed by the boreal and Arctic regions [START_REF] Graven | Enhanced Seasonal Exchange of CO₂ by Northern Ecosystems Since 1960[END_REF][START_REF] Piao | Weakening temperature control on the interannual variations of spring carbon uptake across northern lands[END_REF]. Our findings of a stall in photosynthetic carbon uptake in East Asia and western North America suggest that climate change could have profoundly altered how terrestrial ecosystems over the temperate northern hemisphere responds to external forcing. We might expect a stall in the future for the trends of photosynthesis and SCA at higher northern latitudes (e.g. BRW), if these northern ecosystems become more drought-stressed. Indeed, a weakening positive, or even emerging negative, impact of temperature on SCA has been observed at BRW over the last three decades, and this change in impact of temperature is probably explained by the co-occurring drought with heatwave [START_REF] Peñuelas | Shifting from a fertilization-dominated to a warming-dominated period[END_REF].

Since increasing SCA at MLO was primarily before the 1980s, while contribution of agricultural intensification to SCA increase at MLO was estimated to mainly occur after the mid-1980s (Zeng et al., 2014, it put into question whether agricultural intensification was indeed the dominant driver of SCA enhancement at MLO over past six decades. Our transport experiments based on NBP estimated from both DGVMs and atmospheric inversions do not support agricultural zones dominating SCA increment (Figure 6; Figure S20), though the estimates are partly limited by the prior land flux used in atmospheric inversions and the incomplete representation of agriculture practices in DGVMs. We also indicate that CO 2 fertilization may not always be the most prominent factor driving the SCA trend. Climate change-induced intensified drought stress in temperate northern ecosystems explains a large part of deceleration of SCA at MLO since the 1980s. For high latitude sites, climate change has exerted strong influence on the change of SCA through many other effects, like the structure change of biosphere [START_REF] Forkel | Enhanced seasonal CO 2 exchange caused by amplified plant productivity in northern ecosystems[END_REF][START_REF] Fung | A Hyperventilating Biosphere[END_REF]

Accepted Article

This article is protected by copyright. All rights reserved 2013), reduced residence time [START_REF] Jeong | Accelerating rates of Arctic carbon cycling revealed by long-term atmospheric CO 2 measurements[END_REF], seasonal compensation [START_REF] Liu | Increased high-latitude photosynthetic carbon gain offset by respiration carbon loss during an anomalous warm winter to spring transition[END_REF] and increased winter respiration [START_REF] Commane | Carbon dioxide sources from Alaska driven by increasing early winter respiration from Arctic tundra[END_REF].

Analysis of the SCA trend and its change over the past six decades are one of a few ways to detect changes in the global carbon cycle from long-term CO 2 records [START_REF] Ballantyne | Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years[END_REF][START_REF] Wang | A two-fold increase of carbon cycle sensitivity to tropical temperature variations[END_REF]. These diagnostics from the measurements, and the models' skill in reproducing them, can improve our understanding of the ecosystem response to climate change, and provide useful benchmarks for DGVMs, as long as the impact of atmospheric transport can be accurately considered. Additional long-term records at different latitudes from denser measurement networks over different continents could improve our capability to detect potential changes in the state of the carbon cycle (e.g., [START_REF] Peñuelas | Shifting from a fertilization-dominated to a warming-dominated period[END_REF], reduce uncertainties in the regions and processes that contribute to the change, and enhance the signal to noise ratio for detecting decadal variations of the carbon cycle. They would allow us to more accurately attribute and project the evolution of the carbon cycle in a warmer world with higher atmospheric CO 2 mole fraction, higher surface temperatures and increasing frequency and intensity of droughts [START_REF] Cai | Increasing frequency of extreme El Niño events due to greenhouse warming[END_REF]IPCC, 2013). Changes in the trends of NBP between the two time periods under the effect of CLIM, CO 2 and LU are shown in panels (f), (i) and (l), respectively. In panels (d)-(l), white scatters indicate that the signs of the NBP trend derived from less than eight models are consistent with the BMA trend.

Accepted Article

For better clarity, the sign was plotted on a 1 degree grid with each cell containing four 0.5 degree pixels. We compared the NBP trend between 1959-1984 and 1985-2016 because a piecewise regression [START_REF] Toms | Piecewise regression: a tool for identifying ecological thresholds[END_REF][START_REF] Wang | Spring temperature change and its implication in the change of vegetation growth in North America from 1982 to 2006[END_REF]

(

  NBP clim ) can also contribute the trend change of SCA at MLO (Figure 4f). To better understand the contribution of regions where the climate change induced reversal of the NBP trend to the trend change of SCA at MLO, we performed another transport experiment: T9. In this case, NBP only varied over pixels with NBP clim trend reversal (NBP trend decreases under the impact of climate change) between 1959-1984 and 1985-2016, with the NBP of other pixels kept constant at their 1959 values throughout the 1960-2016 simulation (Figure 4f). The contribution of NBP clim trend reversal regions to the trend change of SCA at MLO can be derived by comparing Simulation T9 and Simulation T2. To quantify potential contribution from the NBP clim trend reversal regions of South America, we performed another transport experiment T9 S.Am to estimate the contribution of NBP trend in NBP clim trend reversal regions of South America to the SCA trend change at MLO. In this transport experiment, similar with Simulation T9, we allowed NBP of south America where NBP clim trend has decreased since the 1980s to vary, and kept NBP of other land area constant at their 1959 values throughout 1960-2016 simulation period. By comparing T9 S.Am and T2, we obtained the contribution from NBP clim trend reversal regions of South America to SCA trend change at MLO.

  In summary, with an ensemble of dynamic global vegetation models coupled with an atmospheric transport model, we perform transport experiments to hindcast the change of SCA and isolate the effects of major factors and regions driving this change. We found the deceleration of SCA since the mid-1980s largely resulted from the response of land carbon cycle to climate change and from changes in atmospheric circulation. Climate change increased SCA at MLO before the 1980s but decreased it afterwards, which is probably associated with the intensified drought stress since the 1980s over the temperate Northern Hemisphere. The critical role of change in atmospheric circulations highlights the long-dismissed necessity to adequately account for changing circulation patterns in understanding carbon cycle change from atmospheric observations and benchmarking of the DGVMs, whose accuracy in hindcasting carbon cycle change remains to be further improved.

FIGURE 1

 1 FIGURE 1 Anomalies of atmospheric seasonal CO 2 amplitude (SCA) at Mauna Loa (MLO). The observed anomaly in the SCA at MLO (black line) was derived from the Scripps CO 2 Program during 1959-2018. The modeled SCA anomaly for 1959-2016 was estimated from a global atmospheric transport model run (Simulation T1) with net biome productivity, fossil fuel CO 2 emission and ocean-atmosphere CO 2 exchange fluxes. Simulated SCA, driven by the NBP from twelve DGVMs, was averaged using two approaches: arithmetic model averaging (AMA) (blue line) and Bayesian model averaging (BMA) (yellow line), respectively. The shaded areas represent the 1 SD uncertainty across models.

FIGURE 2

 2 FIGURE 2 Distribution of observed seasonal CO 2 amplitude (SCA) trends at Mauna Loa (MLO) before and after the 1980s and of relative change of SCA trend between two periods. In panel (a), the histograms indicate the results of SCA trend merging distributions with different breakpoints (1983-1986). In panel (b), relative change of SCA trend is the ratio of calculated difference between SCA trend after the 1980s and SCA trend before the 1980s to SCA trend before the 1980s. Anomalous years are randomly removed in the trend estimation.

FIGURE 3

 3 FIGURE 3 Matrices of seasonal CO 2 amplitude (SCA) trend at Mauna Loa (MLO). The matrices of the SCA trend at MLO during periods with different combinations of start and end year are shown in panel (a) (observations) and (b) (model ensemble with Bayesian model averaging). Panels (c)-(h) show the trend matrices based on BMA at MLO under the effects of six different factors. In each panel, the magnitudes of the trends are shown in the upper left part and the corresponding P values are shown in the lower right part. The individual effects of change in climate ('CLIM') (c), atmospheric CO 2 ('CO 2 ') (d), land use ('LU') (e), fossil fuel ('FF') (f), ocean-air carbon flux ('Ocean') (g), and atmospheric transport ('Wind') (h), on trends of SCA were estimated from transport simulations by an atmospheric transport model fed with land/ocean carbon fluxes.

FIGURE 4

 4 FIGURE 4 Spatial pattern of the trend in net biome productivity (NBP) of the carbon uptake period (May-September) during 1959-1984 and 1985-2016. Panels (a) and (b) show the mean Accepted Article

  of SCA at MLO shows that the trend change occurs in 1984. The shaded area indicates terrestrial regions where the mean surface flux sensitivity is less than 150 ppm per kg C m -2 h -1 , and ocean. The grey areas over land in panels (d)-(l) indicate the regions where NDVI is lower than 0.1.

FIGURE 5

 5 FIGURE 5 Trend of Palmer Drought Severity Index (PDSI) and simulated soil moisture in the carbon uptake period (May-September) during 1959-1984 and 1985-2016. The soil moisture used here is the average of TRENDY models using BMA. The time span in panel (a) and (b) is from 1959 to 1984, and that in panel (c) and (d) is from 1985 to 2016. Change in the trend of PDSI and soil moisture between 1985-2016 and 1959-1984 is shown in panel (e) and (f), respectively. Shaded area indicates regions where mean surface flux sensitivity is less than 150 ppm per kg C m -2 h -1 and ocean. The grey areas over land indicate the regions where NDVI is lower than 0.1. Lower PDSI indicates drier condition.

FIGURE 6

 6 FIGURE 6 Contribution of agricultural intensification to seasonal CO 2 amplitude (SCA) anomaly (b) and trend (c). Spatial distribution of cropland (yellow) is shown in panel (a). The contribution of cropland to SCA anomaly (red line) was estimated with two long-term inversion models, Jena_s57X and Jena_sEXT from the difference between simulation TI and TIC (panel b). The uncertainty of SCA anomaly from inversion models was estimated as the range of two models (b). Observed SCA trend (black bars) and SCA trend contributed by cropland from inversion Accepted Article

LMDZ transport model was nudged to ECMWF analyzed winds.

Regions where mean surface flux sensitivity is larger than 150 ppm per kg C m -2 h -1 .

Regions where the fraction of cropland is larger than 30%.Accepted Article

ACKNOWLEDGEMENTS K. W. and Y. W. contribute equally to this study. This study was supported by the National Key R&D Program of China (2017YFA0604702), the National Natural Science Foundation of China (41530528) and the 111 Project (B14001). This study was granted access to HPC resources of CCRT under allocation A0050102201 made by GENCI (Grand Équipement National de Calcul Intensif).

Accepted Article

This article is protected by copyright. All rights reserved https://doi.org/10.1073/pnas.1014425108 Wei, Y., Liu, S., Huntzinger, D. N., Michalak, A. M., Viovy, N., Post, W. M., … Shi, X. (2014). The North American Carbon Program Multi-scale Synthesis and Terrestrial Model Intercomparison Project -Part 2: Environmental driver data. Geoscientific Model Development, 7, 2875-2893. https://doi.org/10.5194/gmd-7-2875-2014 Wenzel, S., Cox, P. M., Eyring, V., & Friedlingstein, P. (2016). Projected land photosynthesis constrained by changes in the seasonal cycle of atmospheric CO₂. Nature, 538, 499-501. https://doi.org/10.1038/nature19772 Zeng, N., Zhao, F., Collatz, G. J., Kalnay, E., Salawitch, R. J., West, T. O., & Guanter, L. (2014). Agricultural Green Revolution as a driver of increasing atmospheric CO₂ seasonal amplitude. Nature, 515, 394-397. https://doi.org/10.1038/nature13893 Zhao, F., Zeng, N., Asrar, G., Friedlingstein, P., Ito, A., Jain, A., … Zaehle, S. (2016). Role of CO₂, climate and land use in regulating the seasonal amplitude increase of carbon fluxes in terrestrial ecosystems: A multimodel analysis. Biogeosciences, 13, 5121-5137. https://doi.org/10.5194/bg-13-5121-2016

Accepted Article

This article is protected by copyright. All rights reserved

DATA AVAILABILITY STATEMENT

The DGVM simulations are available upon request from S. Stich (s.a.sitch@exeter.ac.uk) or P.

Friedlingstein (p.friedlingstein@exeter.ac.uk). Ocean model simulations are available upon request from C. Le Quéré (c.lequere@uea.ac.uk). Other data supporting the findings of this study are available within the paper. Any additional information may be obtained from the corresponding author upon reasonable request.

Accepted Article

models (red and orange bars) at Mauna Loa (MLO) for two periods, 1959-1984 and 1985-2014 are shown in panel (c). Change of SCA trend indicates trend in 1985-2014 minus that in 1959-1984. The uncertainty of the SCA trend was estimated as the standard error of the linear regression coefficient, and the uncertainty of change in SCA trend was estimated using bootstrap analysis.

The white areas over land in panel (a) indicate the regions where NDVI is lower than 0.1.

Accepted Article

This article is protected by copyright. All rights reserved

Accepted Article

This article is protected by copyright. All rights reserved