Decadal Variability in Land Carbon Sink Eciency Reveals Apparent Trend Reversal After 2009

Background: The climate mitigation target of limiting the temperature increase below 2 °C above the preindustrial levels requires the efforts from all countries. Tracking the trajectory of the land carbon sink eciency is thus crucial to evaluate the nationally determined contributions (NDCs). Here, we dene the instantaneous land sink eciency as the ratio of natural land carbon sinks to emissions from fossil fuel and land-use and land-cover change with a value of 1 indicating carbon neutrality to track its temporal dynamics in the past decades. Results: Land sink eciency has been decreasing during 1957-1990 because of the increased emissions from fossil fuel. After the effect of the Pinatubo eruption diminished (after 1994), the land sink eciency rstly increased before 2009 and then began to decrease again after 2009. This reversal around 2009 is mostly attributed to changes in land sinks in Latin America in response to climate variations. Conclusions: The decreasing trend of land sink eciency in the recent years reveals greater challenges in climate change mitigation, and impacts of climate on land carbon sinks need to be accurately quantied to assess the implementation of climate mitigation policies. change; ENSO:El-Niño-Southern Oscillation; CAMS:Copernicus Atmosphere Monitoring Service; R sc :slope change values; H&N:Houghton and Nassikas; DGVMs:dynamic global vegetation models; MEI:Multivariate ENSO Index; PDO:Pacic Decadal Oscillation; ESA CCI:European Space Agency Climate Change Initiative


Background
The Paris Agreement, aiming at limiting the temperature increase below 2 °C above pre-industrial levels, also aims at a greenhouse gas (GHG) balance between anthropogenic emissions and sinks by the middle of the century. Although the Paris agreement focuses on anthropogenic ux, in reality it is hard to separate the anthropogenic contribution to global sinks [1]. For example, countries have a discretionary option to declare parts of their territory as being under management or not [2], when evaluating or setting the mitigation targets e.g. in the national determined contributions (NDCs). Despite various possible interpretations of what the GHG balance in the Paris Agreement is [1], understanding trends and variations in the global balance between carbon sources and sinks is essential for the evaluation of the NDCs to climate mitigation. The instantaneous land carbon sink e ciency (E) is de ned here as the ratio of 'natural' carbon ux excluding land-use disturbed areas (B) to emissions from fossil fuel (F) and landuse and land-cover change (L), i.e. E = B/(F + L). E = 1 indicates a carbon neutral region where natural sinks offset land use emissions and fossil fuel emissions (a positive sign of carbon uxes being adopted for carbon emissions to the atmosphere for F and L, and for carbon uptake from the atmosphere for B). E integrates trends from emissions and sinks and thus is relevant for assessing regional trajectories with respect to carbon neutrality.
Changes in individual carbon uxes contribute to the variations of E. For example, global annual F emissions have been increasing from an average of 3.04 ± 0.41 Pg C yr − 1 in the 1960s to 9.5 ± 0.4 Pg C yr − 1 in the past decade (2009 − 2018), largely driven by an increase in China [3,4]. Global annual L emissions ranged from 1.0 to 1.8 Pg C yr − 1 since 1959, mainly contributed by tropical regions (South and Southeast Asia, Latin America and Sub-Saharan Africa) [3,5]. However, emissions from L that are already reported under the United Nations Framework Convention on Climate Change (UNFCCC) and the Global Stocktake are generally much lower than L based on the scienti c de nitions of Global Carbon Project (GCP), because they incompletely account for land degradation emissions, do not account for changes in cropland and grassland management intensity, ignore some conversions of carbon rich biomes like tropical forests becoming plantations (e.g. oil palm, rubber). Although the trend of B from 1980 to 2014 is not signi cant in the atmospheric global carbon budget, the net land ux (BL = B + L) increased signi cantly at a rate of 0.037 ± 0.017 Pg C yr − 2 [6].
The land sink e ciency (i.e. E) is in uenced by anthropogenic activities (e.g. fossil fuel emissions, deforestation, afforestation and reforestation) and climate change, including warming trends and variations of El-Niño-Southern Oscillation (ENSO), which dominate annual anomalies of land sinks.
Volcanic eruption and large-scale re events are also important components that regulate E, both locally and globally [7][8][9][10][11]. The rising atmospheric CO 2 concentrations can enhance plant photosynthesis and thus may increase land sinks [12][13][14], although the CO 2 response of carbon sequestration in mature forests may be insigni cant [15]. While deforestation results in carbon emissions and reduces land sinks, secondary regrowth of forests has the potential to enhance land sinks [16][17][18]. Climate conditions like temperature and precipitation have multiple effects on the land carbon uptake. Although global warming has extended the growing season length and thus enhanced vegetation productivity in the northern temperate and boreal regions [19], autumn warming has also led to offsetting carbon losses from northern ecosystems due to respiration increase [20]. ENSO-induced temperature and precipitation variations also strongly impact the carbon cycles in the tropics and play a dominant role in the variability of land sinks [21].
The objective of this study is to characterize the global and regional trajectories of the land carbon sink e ciencies in the past decades using multiple datasets. We identify critical regions in the global trends of E and analyze the potential driving factors. Our analyses focus on the period covered by atmospheric inversion estimates, especially after 1990s when expansion of the atmospheric station network allowed for latitudinal resolution of surface uxes.
These inversions combine xed fossil fuel emission priors, and adjust land and ocean uxes with their prescribed uncertainties to atmospheric CO 2 observations with a transport model to estimate BL and net ocean sinks distributions.

Page 4/17
The temporal coverage and spatial resolution of CAMS (v18.3) are 1979-2018 and 1.875° latitude × 3.75° longitude. The number of stations used in the CAMS increased over time as they became available. The Jena CarboScope inversion (v4.3) used different station networks with ve products, including a network xed to the 1976 con guration (hereafter Jena_s76, 9 stations), the 1981 (Jena_s81, 14 stations), the 1985 (Jena_s85, 21 stations) or the 1993 ones (Jena_s93, 35 stations), respectively. Moreover, we used the Jena CarboScope run sEXTocNEET_v4.3 (hereafter Jena_s57T) using a growing network of 89 stations starting from 1957, but with year-independent degrees of freedom regressing interannual BL variations against variations in air temperature. The spatial resolution is 3.75° latitude × 5°l ongitude. Because prior xed fossil fuel emissions (F) are different in these two inversions ( Figure S1b), directly impacting posterior BL estimates, the BL data from the Jena inversions were adjusted to a common F value, as done in Peylin et al. [2013] and Thompson et al. [2016]. For global F, we used the values from global carbon budget [Global Carbon Project, GCP, 27] and for regional F, we used the values from CAMS which is closer to those assessed by the global carbon budget. All bunker fuels being considered as a surface fossil CO 2 source distributed proportionally to national emissions shares of the global total [4,27].
We also used estimates of B aggregated at the scale of both hemispheres over 1994-2013, using a twobox inversion and data from the two longest CO 2 monitoring stations from South Pole and Mauna Loa, with ocean sinks from an ensemble of ocean biogeochemical models [28].

Natural land sink as residual of GCP's annual global carbon budget
We calculated annual global B over 1959-2018 from the global carbon budget [26] as a residual of F, L from bookkeeping models, global atmospheric CO 2 growth, and ocean sinks from an ensemble of ocean biogeochemical models as well. Deriving B as a residual term has recently been replaced by explicit simulations with dynamic global vegetation models (DGVMs) in recent versions of the annual global carbon budget [29]. Nevertheless, we here stay with the residual approach because the land sink from DGVMs added to ocean sinks and F emissions do not match the CO 2 growth rate, with an imbalance ranging from − 1.75 to 1.96 PgC yr − 1 during 1959-2018, so that B from this approach is not consistent with atmospheric data [27,30]. We considered B from the TRENDY-V8 DGVMs used in Friedlingstein et al.
[2019] as a sensitivity test and acknowledge that potential error terms in the other uxes will be attributed to B with residual approach.
Land-use and land-cover change ux from bookkeeping models and DGVMs We mainly used L from the Bookkeeping of Land Use Emissions model [BLUE,31] because this dataset is grid-based and updated to the latest year. The BLUE model used gridded LUC data from the Land Use Harmonization dataset [32]. In addition, we used L from three other sources as sensitivity tests: 1) L from the bookkeeping model by Houghton and Nassikas [2017, hereafter H&N], 2) L from DGVMs in , and 3) L from the OSCAR compact Earth system model that emulates the carbon cycle of TRENDY-V7 DGVMs [OSCAR,33]. L from H&N is estimated based on country-level response curves of carbon pools for different LUC types and FAO/FRA forest data [5]. It should be noted that although H&N and BLUE are used as equally likely in the GCP, H&N is not used here for the prime analysis due to its ending in 2015 and being lack of spatially explicit values after that year (GCP extended the results of H&N to 2018 only on global scale). In the GCP, DGVMs performed two simulations using different settings: S2 with varying CO 2 and climate but time-invariant preindustrial land use maps, and S3 with annually updated CO 2 , climate and land use maps. L is thus the net biome productivity (NBP) difference between S2 and S3, which includes a foregone land sink in S2, leading to a component of L known as "loss of additional sink capacity" (LASC) which does not exist in observation-based estimates and in bookkeeping models [34,35]. OSCAR embeds processes and parameters calibrated using outputs from DGVMs and calculates L using a bookkeeping method [33]. While L from DGVMs include LASC, L from OSCAR does not and can be compared with other bookkeeping models based on observations of carbon densities [33].
All carbon ux datasets used in this study are summarized in Table S1.

Data analysis
Gridded datasets with different spatial resolutions were resampled to 1°×1°. We only focused on the period covered by inversion data (1957-2018) and used annual mean values of the carbon uxes and climate variables. Because of the strong interannual variability ( Figure S1), we applied a 5-yr moving window to carbon uxes to better detect trends. In the atmospheric inversion output, however, it is impossible to distinguish B and L separately. We thus used L simulated by BLUE to calculate B. We also calculated trends of each individual ux on the regional scale for 11 regions (see region division in Figure   S2) to elucidate which region dominates the trend of E. Piecewise regression ("Segmented" package in R) was applied on global and regional trends to detect breakpoints. In addition, a linear least-square regression was used to calculate the trends during 1957-1990 and the trends before and after a detected breakpoint for each individual ux.
Because the global trend in B/(F + L) is not equal to the sum of regional trends, we used a method of removing-one-region at a time to analyze the regional contributions to the breakpoint of global land sink e ciency. Speci cally, if there is a breakpoint detected in the global signal, the slopes before and after the breakpoint are s 1 and s 2 . After removing B, F and L uxes from one region, we assume that B/(F + L) from the sum of the other regions still shows a similar breakpoint, but the slopes before and after the breakpoint change to s 1 ' and s 2 '. We de ne the contribution of a region as a ratio of slope change (R sc ) : A positive value of R sc indicates that this region strengthens the slope reversal, i.e. making the contribution to the global breakpoint more signi cant. Conversely, a negative value indicates the region weakens the slope reversal.  Figure S1). After the Pinatubo eruption period, E increases until around 2009 and then decreases afterward (Fig. 1b). This trend is detected by all datasets, and the detected breakpoints range from 2008 to 2010 (p < 0.01). In the following parts, we focus mainly on E after the eruption of Pinatubo because more datasets with better observation constraints are available during this period (i.e. from 1994-2018). The original ux B and F are generally increasing during 1957-2018 ( Figure S1). B from different datasets is roughly consistent, but L from different dataset shows large variations over 1957-2018 owing to the various methods and input datasets ( Figure S1).

Regional contributions
In general, the global breakpoint around 2009 remains detectable after removing contributions from one individual region in both inversions and is all signi cant (p < 0.1) except when removing Latin America in Jena_s93 ( Figure S3-S8). The regional contributions could be clearly re ected by their R sc values (Fig. 2a). The mean values of R sc show that Latin America contributes most to the slope reversal with the largest positive R sc values (i.e. enhancing effect on the slope reversal) while East Asia shows largest opposite effects with the most negative R sc values. The negative contributions from North America, Europe and Middle-East are also consistent in all datasets (R sc < 0 for all). However, Former Soviet Union shows more positive contributions than Africa in Jena_s81, Jena_s85 and Jena_s93 (Fig. 2a) but opposite contributions are found in CAMS (Fig. 2a). East Asia, North America and Europe all have high F and the most negative R sc values, which proved that the reversal was not due to a change of F in any of these large F emitting regions.
Given the robustness of the global breakpoint after removing each individual region ( Figure S3-S8), we further removed 2 regions each time to examine whether the breakpoint would disappear. We nd that several 2-regions combinations can make the breakpoint insigni cant (p > 0.1) or disappear (Fig. 2b).
Among these combinations, Latin America is the most frequently identi ed region and Former Soviet Union ranks second. However, these two regions are not sensitive regions in CAMS. In CAMS, only removing the combination of Africa and Southeast Asia can change the breakpoint. After removing any combination without Latin America, Africa, Former Soviet Union and Southeast Asia, the breakpoint still exists, con rming the minor or even opposite contributions of these regions.
Although the two-box inversion cannot detect the breakpoint due to the short time series after 2009, it con rms that the increasing trend before 2009 is contributed by the Southern Hemisphere ( Figure S9).

Individual uxes in each region
Because it is partly ambiguous to simply decompose the trends in global E = B/(F + L) into trends in regional E, we calculate the global and regional trends in each individual ux before and after 2009 from 1996-2016 in different datasets (Fig. 3). Globally, F and L are increasing before and after 2009 while B is increasing only before 2009. After 2009, the growing trend of F weakens while the growing trend of L increased, albeit with considerable variability among estimates ( Figure S1c). The trend in B changes from positive to negative after 2009. Therefore, the reversal of trends in E are driven predominantly by the trend in B and, to a lesser extent, by the trend in L. The trend in F, on the other hand, shows a more slowly increasing rate after 2009 than before 2009, making the E value more positive after 2009. Thus, the trend in E is primarily driven by land sink variability (e.g. B), and F emissions have an opposite effect but with small magnitude.
The global trend in each carbon ux is the sum of regional trends. The global trend in F is mostly driven by trends in East Asia where the growth of fossil fuel emissions slows down after 2009 compared to that before 2009 (Fig. 3). A slow-down in F growth also appears in Africa, Middle East, Oceania and Southeast Asia with a smaller magnitude. In North America, however, F shows a small magnitude of decreasing

Discussion
We discovered a decreasing trend of E from 1957-1990 which is in line with the previous nding that the proportion of carbon emissions remaining in the atmosphere (airborne fraction) was increasing during this period [42]. After the Pinatubo eruption period (after 1994), E increased again and then decreased after 2009. To verify the robustness of this breakpoint, we did the following sensitivity tests: 1) using L from other three data sources (H&N, DGVMs and OSCAR), B from DGVMs and B from the residual of different ocean ux estimates; 2) calculating E using all the combinations (28 in total) of B from seven estimates and L from four estimates; 3) removing adjustment of fossil fuel emissions on Jena data; 4) applying different moving average methods; and 5) masking the latest El Niño event.
In general, our breakpoint is robust regardless of different choices of datasets and methods. Although there are large differences in L among different datasets ( Figure S1c), the absolute values of L are small, and the breakpoint detection in E is robust regardless of different L choices (p < 0.01, Figure S11). This breakpoint is, however, not re ected using B from DGVMs, mainly due to the large interannual variability of B before 2009 ( Figure S12). On the other hand, the decreasing trend of E after 2009 in DGVMs ( Figure  S12) is consistent with the trend detected by inversion data. A previous study also found the global difference in land sink between DGVMs and inversion datasets agreed well with the budget imbalance in the global carbon budget [30], indicating a possible bias in the land sink simulated by DGVMs. We also test the ocean sink from each individual ocean model estimate and pCO 2 -based product of the global ocean sink reported in the global carbon budget, instead of the ensemble mean, to calculate the residual B in the global carbon budget, and the breakpoint of global E still exists (p < 0.01, Figure S13). Using the mean E from 28 combinations of B and L estimates, the detection of breakpoint remains signi cant (p < 0.01, Figure S14). With more observation data constraint, the 1-σ uncertainty range of E is smaller after the Pinatubo eruption ( Figure S14). Without the adjustment of fossil fuel emissions on Jena data (see Sect. 2.1), the global pattern remains consistent (p < 0.01, Figure S15). To evaluate the in uence of moving average methods, we applied 3-yr moving average on the time series, and signi cant global breakpoint around 2009 in E is found in each dataset (p < 0.1, Table S2). Even using the original annual values without moving averages, the breakpoint around 2009 is still detectable although not signi cant (Table S2). After replacing the B in the strong El Niño years (2015 and 2016) with the averaged B of 2014 and 2017, the breakpoint detection is still signi cant after 5-yr moving average in CAMS and 3 Jena datasets (p < 0.1, Table S2), indicating that this reversal is not completely caused by the latest strong El Niño event during 2015-2016. These sensitivity tests further con rm the robustness of the breakpoint and support that the smaller increasing trend of B is the dominant factor determining the reversed trend in land sink e ciency after 2009.
Compared to the period before 2009, the acceleration of L growth in Africa and East Asia together with the weakening or relatively stable trend in B in the tropics, especially in Latin America, result in the weakening of E after 2009 (Fig. 3). This is consistent with the dominant role of tropical regions in the interannual variability of the global carbon cycle [43][44][45][46].
Climate variations play an important role in the reversed trend in land sink e ciency after 2009. In fact, 5yr moving averages of MEI and PDO index both show a strong increasing trend since 2009 ( Figure S16), which in uences the land carbon sink. However, the impacts of ENSO and PDO on E seem to be small ( Fig. 1) during the increasing phase of ENSO and PDO in the previous cycle (2000)(2001)(2002)(2003)(2004), probably because of reduced intensity and shorter duration ( Figure S16). In fact, if we apply multiple linear regression using annual precipitation, temperature and CO 2 concentration to predict B in tropics and use the predicted B to calculate tropical E, the same breakpoint in the predicted E trend around 2009 is still signi cant ( Figure S17a), indicating these factors can largely explain the trend reversal. Note that MEI is strongly correlated with temperature and precipitation in the tropics ( Figure S17b). Previous studies also found that climate variations caused by e.g. El Niño may in uence the tropical carbon cycle through different processes in different regions [21].
Climate impacts on land sink may need to be accurately quanti ed in the future. It has been proved that extreme El Niño events that strongly reduce tropical land carbon sinks are expected to be more frequent due to future greenhouse warming [47], but the trend of El Niño/La Niña intensity still remains unclear. While El Niño/La Niña cycles affect E, El Niño impacts are probably compensated by recovery uxes in the subsequent years [48]. If they can fully offset each other, ENSO impacts on E may be negligible in the long term. On the other hand, if there is an anthropogenic ngerprint in trends in El Niño/La Niña intensity [49], climate change impact on natural land sink need to be considered to de ne mitigation goals compatible with the Paris agreement. In addition, PDO and other low frequency variability patterns might affect climate at the time scales these mitigation policies should be acting (the next few decades). They could either result in additional CO 2 in the atmosphere (e.g. by imposing more drought/higher temperatures) amplifying the impacts or reduce it temporarily (if they would lead to some cooling and wetter conditions). Their impacts on the natural sink thus need to be accurately quanti ed to avoid a false sense of implementation progress when assessing climate mitigation policies. Currently, there is not enough evidence to identify the most likely of these two possibilities. Nevertheless, it is necessary to track the e ciency of natural sink dynamics.
Saturation of land carbon sinks could also contribute to the reversal of trends in land sink e ciency after 2009. Processes that regulate land sink driven by atmospheric composition change (e.g. CO 2 fertilization, nitrogen deposition), climate change (e.g. rising temperature) and LUC (e.g. forest regrowth, woody encroachment) are unlikely to be sustained permanently. Land carbon sinks are thus likely to decrease as the terrestrial carbon storage saturates. The saturation of the land sink is already indicated by a network of Amazonian forest plots [46] and modelling studies [50]. However, some of this potential saturation may be offset by secondary forest regrowth in the tropics [16,51]. Forest area gain in the tropics during 2009-2013 is lower than over the period before 2009 but became higher in the recent period after 2013 ( Figure S18). Linking forest area gain to the net land sink, however, remains challenging due to the uncertainties in the forest gain detection from remote sensing, biomass growth and soil carbon dynamics. Therefore, the contribution of legacy land sink from forest gain needs to be further investigated with emerging evidence.

Conclusions
Our study de ned the land sink e ciency and explored its trends from 1957-2018 using multiple datasets. We found that before the Pinatubo eruption, the land sink e ciency was generally decreasing due to the increasing in fossil fuel emissions. After the Pinatubo eruption, a trend reversal around 2009 was observed because of the changed land sink in Latin America. Our results highlight the importance of taking into account variations of the natural land sink, which is not explicitly included in the Paris agreement. The decreasing trend of land sink e ciency in the recent years reveals greater challenges in the mitigation of climate change and reinforces the need for more efforts to reduce anthropogenic emissions from fossil fuel burning and LUC. Global Carbon Budget data used in this study can be accessed from https://doi.org/10.18160/gcp-2019 [27]. Jena CarboScope inversion data is available at http://www.bgc-jena.mpg.de/CarboScope/.

Competing interests
The authors declare that they have no competing interests.