Evidence of Carbon Uptake Associated with Vegetation Greening Trends in Eastern China

Abstract

Persistent and widespread increase of vegetation cover, identified as greening, has been observed in areas of the planet over late 20th century and early 21st century by satellite-derived vegetation indices. It is difficult to verify whether these regions are net carbon sinks or sources by studying vegetation indices alone. In this study, we investigate greening trends in Eastern China (EC) and corresponding trends in atmospheric CO₂ concentrations. We used multiple vegetation indices including NDVI and EVI to characterize changes in vegetation activity over EC from 2003 to 2016. Gap-filled time series of column-averaged CO₂ dry air mole fraction (XCO₂) from January 2003 to May 2016, based on observations from SCIAMACHY, GOSAT, and OCO-2 satellites, were used to calculate XCO₂ changes during growing season for 13 years. We derived a relationship between XCO₂ and surface net CO₂ fluxes from two inversion model simulations, CarbonTracker and Monitoring Atmospheric Composition and Climate (MACC), and used those relationships to estimate the biospheric CO₂ flux enhancement based on satellite observed XCO₂ changes. We observed significant growing period (GP) greening trends in NDVI and EVI related to cropland intensification and forest growth in the region. After removing the influence of large urban center CO₂ emissions, we estimated an enhanced XCO₂ drawdown during the GP of −0.070 to −0.084 ppm yr⁻¹. Increased carbon uptake during the GP was estimated to be 28.41 to 46.04 Tg C, mainly from land management, which could offset about 2–3% of EC’s annual fossil fuel emissions. These results show the potential of using multi-satellite observed XCO₂ to estimate carbon fluxes from the regional biosphere, which could be used to verify natural sinks included as national contributions of greenhouse gas emissions reduction in international climate change agreements like the UNFCC Paris Accord.

1. Introduction

Greening trends appear over much of the global land area [1], indicating changes in ecosystem function at regional scales. Atmospheric carbon dioxide (CO₂) fertilization, temperature, and precipitation changes account for most of the greening of the Northern Hemisphere (NH) over the past two decades [2,3,4]. Terrestrial greening in cold climate zones is considered primarily due to warming temperature [5]. Greening may also be caused by the precipitation changes in cool-arid systems like Northern China and Mongolia [6]. In addition to changing climate, vegetation greening can also be related to forest regrowth following the abandonment of agricultural lands or plantations [7]. The drivers of greening are fairly well understood and correspond to changes in photosynthetic CO₂ uptake [2,4]. However, relating these changes in vegetation indices to changes in net carbon fluxes is challenging, because they are largely insensitive to CO₂ release by respiration fluxes [8].

Large-scale land management in certain areas, including agricultural expansion and intensification and afforestation and reforestation efforts, is expected to alter net carbon uptake at the local to regional scales. China’s land management policies have promoted both of these practices in recent decades [9,10,11,12,13]. Globally, agricultural production has intensified from the last century to maintain food production for population growth [14,15]. Mixed influences on carbon fluxes have been shown, for example, increased carbon uptake due to longer growing seasons [16] or intensification [17] but stronger carbon emissions from the conversion of natural systems to cultivated crops [18]. Sustainable cropland intensification is necessary for China using available arable land and water [10]. If done properly, it has the potential to reduce greenhouse gas emissions, through careful management of soil, nutrients, and water [19].

Forests are significant contributors to the global terrestrial carbon sink [20], absorbing approximately one-quarter to one-third of the anthropogenic CO₂ emissions. Planted forests (afforestation and reforestation) may be used to offset anthropogenic CO₂ emissions [21]. Planted forests across China increased significantly at the end of 20th century, increasing the carbon density from 15.3 to 31.1 mega-grams per hectare [22]. Forests planted by China’s Grain for Green Program, started in 1999, are estimated to have increased natural carbon sequestration, thereby mitigating CO₂ emissions from fossil fuel and cement production [12,13] by about 3–5% of China’s annual carbon emissions (relative to 2010 emissions), which may continue for next 40 years as these forests continue to grow [11]. Even though the precise forest area change is difficult to estimate, the terrestrial carbon sink increase over China is generally accepted [9,23], with a net sink in the range of 0.26 Pg C per year absorbing 28–37% of the country’s cumulated fossil carbon emissions during the 1980s and 1990s [9,23,24].

In addition to vegetation indices, atmospheric CO₂ concentrations using column-averaged CO₂ dry air mole fraction (XCO₂) provide an alternative way to investigate terrestrial biosphere carbon flux changes [25]. XCO₂ represents the total amount of CO₂ in a vertical column through the atmosphere, and therefore is less sensitive to changes in boundary layer height and surface winds that impact a concentration measurement near the surface without necessarily a corresponding change in local CO₂ fluxes. Regional CO₂ fluxes have been inferred using satellite observed XCO₂ and atmospheric transport models, even though some challenges exist, like satellite observations coverage, systematic errors of the retrievals, and systematic errors of the transport models [26,27]. XCO₂ seasonal cycle changes are strongly correlated with vegetation indices from satellite observations from 2009 to 2014 [28,29]. Extremely high XCO₂ anomalies are usually accompanied by events such as extreme droughts or floods that reduce vegetation productivity [30,31]. Parazoo et al. [32] found a decrease in tropical XCO₂ corresponding to the gross primary production (GPP) increase in the wet season and a XCO₂ increase associated with strong biomass burning in dry season in Southern Amazonia. These studies show that XCO₂ observations started in 2002 are just beginning to yield insight into variability and trends in surface CO₂ fluxes as the records now extend to just over 15 years.

Regional CO₂ fluxes have been inferred from satellite-observed XCO₂ and in-situ surface observations using atmospheric transport models [26,27,33]. Incorporating XCO₂ into atmospheric inversions improves upon sparse and location-biased in-situ surface observations and yields estimates that are less sensitive to atmospheric transport model uncertainty [34,35,36]. It has also been noted that including XCO₂ constraints in atmospheric inversions shifts regional biospheric uptake from the tropics to the boreal latitudes [33]. Surface flux estimates from atmospheric inversions are most robust across large scales, like the Transcom regions [37] and diverge at sub-continental scales and smaller [35]. This is due largely to uncertainty in atmospheric transport models and the typically sparse surface CO₂ constraints, which may attribute observed atmospheric CO₂ anomalies to many different source/sink spatial distributions. In this analysis, we assume that observed localized spatial XCO₂ anomalies are likely more sensitive to local surface CO₂ fluxes in that specific region than flux estimates derived from atmospheric inversions. This is because there is finer spatial resolution in XCO₂ observations than estimated inversion CO₂ fluxes. Therefore, sub-continental-scale CO₂ fluxes may be inferred from XCO₂ anomalies directly.

Eastern China (EC) is one area with high densities of vegetation cover and greening trends. Even though there are some modeled carbon sequestration estimates of China’s Grain-to-Green Program [11,12,13], the magnitude of CO₂ uptake has not been validated at the regional scale. This study aims to quantify recent trends in regional carbon uptake by quantifying changes in satellite-observed XCO₂ during the growing period (GP). We adopt various satellite-observed vegetation parameters to characterize the seasonality and locations of increased vegetation photosynthetic activity in EC and examine XCO₂ trends in those regions specifically.

2. Data and Methods

Eastern China (EC) mainland, a large carbon storage area in the mid-latitudes [24], was selected for this case study evaluation of carbon uptake trends because it shows persistent vegetation greening (Section 3.1) (outlined in red in Figure 1). Land cover consists of forest, irrigated and rain-fed cropland, and grassland in south (<34°N) and north (>=34°N) parts of EC.

Datasets used in this research are summarized in

Table 1. Datasets used in this study.

Name	Sources	Time Res.	Space Res.	Period	Reference
NDVI/EVI	MOD13C1 v6	16 days	0.05 deg.	2003.01–2016.12	[38]
NDVI	GIMMS 3g	15 days	1/12 deg.		[39]
LAI	MCD15A3H v6	4 days	500 m		[40]
LAI	NOAA CDR v4	1 day	0.05 deg.		[41]
GPP	MOD17A2H v6	8 days	500 m		[42]
Global maps of XCO2	SCIAMACHY BESD-XCO2 GOSAT ACOS-XCO2 OCO-2 XCO2	8 days	1 deg.	2003.01–2009.05 2009.06–2014.08 2014.09–2016.05	[43,44] [45] [46]
XCO2/CO2 flux	CarbonTracker 2017	3 h	2 × 3/1 × 1 deg.	2003.01–2016.12	[47]
XCO2/CO2 flux	MACC-III	3 h	1.875 × 3.75 deg.		[48]
XCO2	GEOS-Chem v11.4	3 h	2 × 2.5/0.5 × 0.625 deg.		[49]
TCCP	AVHRR VCF	Yearly	0.05 deg	2003.01–2016.12	[50]
Fossil fuel CO2 Emission	ODIAC 2017	monthly	1 deg.		[51]
LST	MCD11C2 v6	8 days	0.05 deg.		[52]
Precipitation	TRMM 3B42	3 h	0.25 deg.		[53]

. They include four categories: (1) Vegetation parameters, including normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), leaf area index (LAI), and gross primary production (GPP); (2) XCO₂ time series based on observations from multiple greenhouse gas satellites, including Scanning Imaging Absorption spectrometer for Atmospheric Chartography (SCIAMACHY) onboard the Environmental Satellite (ENVISAT), Greenhouse Gases Observing Satellite (GOSAT), and Orbiting Carbon Observatory 2 (OCO-2); (3) CO₂ fluxes and corresponding modeled XCO₂ concentrations from two inversion models, CarbonTracker (CT) 2017 and Monitoring Atmospheric Composition and Climate (MACC-III), and one forward model GEOS-Chem v11.4; (4) supporting independent datasets, including tree canopy cover percentage (TCCP) and climate factors like land surface temperature (LST) and precipitation, are described in more detail in Section 2.4.

2.1. Trends in Regional Vegetation Indices

2.1.1. Multiple Vegetation Parameters

Satellite observed vegetation indices, included normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) from Moderate-Resolution Imaging Spectroradiometer (MODIS) (MOD13C1 v6) [38] and NDVI from Advanced Very-High-Resolution Radiometer (AVHRR) (GIMMS 3g) [39] were used for greening trend analysis. Leaf area index (LAI), indicating vegetation canopy, from MODIS (MCD15A3H v6) [40] and AVHRR (NOAA CDR v4) [41], were also used for estimating vegetation growth. In addition, GPP (MOD17A2H v6) [42] from MODIS was also used for identifying the trend of photosynthetic productivity.

MOD13C1 v6 was downloaded from National Aeronautics and Space Administration (NASA: https://search.earthdata.nasa.gov/). GIMMS 3g was from National Oceanic and Atmospheric Administration (NOAA: https://data.nodc.noaa.gov/). MCD15A3H v6, MOD17A2H v6, and NOAA CDR v4 were retrieved from Google Earth Engine [54]. NDVI/EVI from MOD13C1 v6, LAI, and GPP were resampled to temporal and spatial resolution of 16 days and 1 × 1 degree using bilinear interpolation. NDVI from GIMMS 3g was also integrated to 1 × 1 degree.

2.1.2. Growing Period Identification

Vegetation indices are widely used to characterize phenology [55,56,57]. We adopted an exponential and logistical smoothing function [57] to fit EVI time series to quantify the beginning and end of strong vegetation activity as the growing period (GP), shown in Equation (1).

$$y\left(t\right)=\frac{c}{1+e^{a+bt}}+d$$

(1)

where t is time in days, $y\left(t\right)$ is the EVI value at time t, a and b are fitting parameters, and c and d were used to constrain the range of EVI values. The spatial/temporal resolution of EVI used here is in 1 × 1 degree/16-day. The start/end of GP for each grid cell can be defined as minimum/maximum value of derivation of the curvature of the function fitted to either the early season data or the late season data separately (Figure 2, the period between two red dashed lines). We excluded grid cells with poorly constrained GP from our analysis including those starting after July, ending before April, or extending for longer than 12 months as these are unrealistic results. The start, end, and length of the GP in each grid cell used for the trend analysis shown in Figure S1. The GP length retrieved in EC was 149 ± 36 days, which started from day of year 126 ± 27 and went to day of year 275 ± 16 with longer GP lengths in the south and shorter in the north (Figure S1).

2.1.3. Vegetation Indices Trend Calculation

We analyzed the trends of multiple vegetation parameters, including NDVI, EVI, and LAI, with a spatial resolution of 1 × 1 degree and time resolution of 15–16 days. After the GP was uniquely defined for each grid cell (Section 2.1.2), the mean of the vegetation indices over the fixed GP length for each year in each grid cell was calculated. Interannual trends in mean GP vegetation indices were defined using least-square linear fitting, and the Mann–Kendal test was used to determine statistical significance of the linear trends in each grid cell.

2.2. Trends in Regional XCO₂ Summer Drawdown

2.2.1. Satellite-Observed XCO₂

Multiple satellite-observed XCO₂ products derived from SCIAMACHY, GOSAT, and OCO-2 platforms were integrated to obtain the longest and most complete record possible [58]. SCIAMACHY on-board the European ENVISAT was launched in 2002, with a temporal resolution of 36 days and spatial resolution of 60 km, available from 2003 to 2012 [59]. GOSAT was launched in 2009 by Japan, with a temporal resolution of 3 days and spatial resolution of 10 km, available from 2009 to the present [60]. NASA’s OCO-2 was launched in 2014, with a temporal resolution of 16 days and spatial resolution of 1.29 × 2.25 km at nadir, available from 2014 to present [61]. SCIAMACHY XCO₂ was retrieved by the full-physics based Bremen Optimal Estimation–DOAS (BESD) algorithm (v02.01.01), referred to as BESD-XCO₂ [43,44] and was downloaded from the European Space Agency (ESA, http://www.esa-ghg-cci.org/). GOSAT was retrieved by the Atmospheric CO₂ Observations from Space (ACOS) project team (v7.3), referred to as ACOS-XCO₂ [45]. OCO-2 was retrieved by the OCO-2 team (v7r), referred as OCO-XCO₂ [46]. ACOS-XCO₂ and OCO-XCO₂ were downloaded from NASA (https://disc.gsfc.nasa.gov/). XCO₂ from these different observations during 2003 to 2016 was integrated and mapped, described in Section 2.2.2, into global maps of XCO₂ (referred as GM-XCO₂) with temporal and spatial resolution of 8 days and 1 × 1 degree.

2.2.2. Merging and Gap-Filling XCO₂ Data Products (Integration and Mapping)

In order to construct a long-term XCO₂ dataset, BESD-XCO₂ from January 2003 to May 2009, ACOS-XCO₂ from June 2009 to August 2014, and OCO-XCO₂ from Sep 2014 to March 2016 were used for XCO₂ integration and mapping. XCO₂ integration and mapping is the extension of the spatio-temporal geostatistical mapping method applied to multiple satellite observations. Corrections were made to account for differences in vertical measurement sensitivity, differences in satellite overpass times and fields of view, and observation gaps. Steps included: (1) Applying an averaging kernel correction based on vertical profiles from CarbonTracker, (2) normalizing overpass times and unifying spatio-temporal scales of multi-satellite XCO₂ data, and (3) global land mapping of XCO₂ based on spatio-temporal geostatistics.

In the first step, we corrected for differences in altitude sensitivity of the sensors and normalized the measurements using averaging kernels and a-priori vertical profiles of CO₂ from CarbonTracker (CT-XCO₂) selected as the common proxy for adjustment (Equation (2)) [43,62].

$$XC{O}_2{}_{adj,t}=XC{O}_2{}_{ret,t}+h^T\left(I-A\right)\left(X_{M,t}-X_{a,t}\right)$$

(2)

where $XC{O}_2{}_{adj,t}$ is the adjusted XCO₂ value at observation time t, $XC{O}_2{}_{ret,t}$ corresponds to retrieved XCO₂ of BESD, ACOS, or OCO-2, h and A are the pressure weighting functions and column-averaging kernels from the retrieval algorithm, respectively, I is an identity matrix, $X_{M,t}$ is the common priori CO₂ profile from CT, and $X_{a,t}$ is a priori CO₂ profile from satellite retrievals.

Considering the diurnal variability in XCO₂ and differences in satellite overpass times, step 2 unified XCO₂ from multiple satellites to a common observing local time and spatial resolution. We introduced a scale factor derived from the diurnal variation of XCO₂ using CT-XCO₂ to normalize the overpass time of XCO₂ retrievals from these three satellites to a reference time (Equation (3)).

$$XC{O}_2{}_{nor,t}=\frac{XC{O}_2{}_{CT,rt}}{XC{O}_2{}_{CT,t}}\ast XC{O}_2{}_{adj,t}$$

(3)

where $XC{O}_2{}_{nor,t}$ is the normalized satellite-observed XCO₂ at the reference time (rt: 13:00, local time); $XC{O}_2{}_{adj,t}$ is the adjusted satellite XCO₂ derived from Equation (2) at satellite overpass time t; $XC{O}_2{}_{CT,rt}$ and $XC{O}_2{}_{CT,t}$ are the XCO₂ from CT model at the reference time rt and satellite overpass time t, respectively.

Additionally, we unified the different field of views (30 × 60 km, diameter of 10.5 km and 2.25 × 1.5 km, respectively) and repeat cycles (36 days, 3 days, and 16 days, respectively) of SCIAMACHY, GOSAT, and OCO-2 before using a gap-filling method. In this study, we used a uniform spatio-temporal scale of 30 × 30 km as a space unit and 8-day as a time unit to average $XC{O}_2{}_{nor,t}$ of BESD-XCO₂, ACOS-XCO₂, and OCO-XCO₂ respectively. This uniform spatio-temporal scale was introduced to aggregate different XCO₂ data products in order to maintain as much as possible the spatial (regional and national) and temporal (monthly and seasonally) variation of XCO₂ observed by these three satellites. The integrated XCO₂ dataset in 30 × 30 km and 8 days was generated from 2003 to 2016, hereafter referred to as integrated-XCO₂. The data quantity and trends of integrated-XCO₂ from 2003 to 2016 is shown in Figure S2.

In step 3, we applied a gap-filling interpolation method to the integrated XCO₂ to produce continuous global XCO₂ maps in time and space. The spatio-temporal geostatistical mapping method, GM-XCO₂, applied in GOSAT observations has been previously used in global carbon research studies [28,31,63,64]. Considering the spatial connectivity used in variogram modelling, we divided global land area (within 45°S–60°N) into five continents. As different atmospheric transport exists between hemispheres, Africa was separated into North Africa and South Africa. As a result, modelling the spatio-temporal correlation structure using variograms of XCO₂ in the gap-filling method was implemented for each of the six land regions: North America, South America, North Africa, South Africa, Europe, and Asia and Oceania, as shown in Figure S3. Detailed information about the spatio-temporal interpolation method can be found in Zeng et al. [63,64]. GM-XCO₂ from Europe and Asia was used in the XCO₂ trend analysis over EC.

GM-XCO₂ data was evaluated using cross-validation (Figure S4) and by comparing to TCCON measurements (Figure S5) [28,31,64]. The predicted XCO₂ (GM-XCO₂) shows a significant linear relationship with the integrated-XCO₂ with R² = 0.97 (p-value < 0.01). Good agreement of integrated-XCO₂ and the gap-filled GM-XCO₂ is shown by the mean predicted error, mean absolute predicted error, and standard deviation of the error equal to −0.02, 0.99, and 1.38 ppm, respectively.

2.2.3. Model Simulations of XCO₂ and CO₂ Fluxes

Estimated surface CO₂ fluxes and XCO₂ concentrations from the GEOS-Chem, CarbonTracker, and MACC models were used for two purposes: (1) To create estimates of XCO₂ background increase in the study area (Section 2.2.4), and (2) estimate relationships between local surface CO₂ fluxes and XCO₂ variability (Section 2.3). CarbonTracker (CT) 2017 is a global inversion model from NOAA providing atmospheric CO₂ mixing ratio distributions and inferred surface CO₂ fluxes with a time resolution of 3 h [47]. It uses the Carnegie–Ames–Stanford Approach (CASA) land surface model to estimate prior flux estimates and Transport Model 5 (TM5) to represent atmospheric circulation [65]. Optimized posterior surface flux components were solved for with a spatial resolution of 1 × 1 degrees including fossil fuel, fire, ocean, and biosphere. Each surface flux was tagged to create associated CO₂ concentrations in 25 vertical atmospheric levels with spatial resolution of 2 degrees in latitude and 3 degrees in longitude. The atmospheric pressure of the different levels used was from the CT CO₂ total files. Contributions to XCO₂ from each component and background concentration were calculated from the vertical profiles based on the pressure-averaged method [66]. MACC-III is also an inversion model from the European Space Agency (ESA) providing global XCO₂ and net CO₂ fluxes with a temporal and spatial resolution of 3 h and 1.875 degrees latitude and 3.75 degrees longitude [48]. Both inversion model results have been compared to the FLUXCOM product in Asia [67], but as both are prone to different types of errors and uncertainties, there is no way to provide a robust validation of these model products.

In order to separate local biospheric CO₂ flux influence in Eastern China from remote biospheric and other non-biospheric CO₂ flux influences on that region, we created estimates of various CO₂ flux components using transport model simulations. These simulation results were later used to remove other non-local and non-biospheric influences from observed GM-XCO2 trends using detrending methods. GEOS-Chem is a global 3-D transport model for atmospheric trace gas simulation driven by meteorological and surface flux input [49]. We adopted GEOS-Chem version 11.4 driven by MERRA-2. In the GEOS-Chem global CO₂ simulation, we used the CarbonTracker fossil fuel, ocean exchange, and biosphere posterior CO₂ fluxes to replace the standard GEOS-Chem surface CO₂ flux input. The biomass fire emissions were taken from the Global Fire Emission Database version 4.1s (GFEDv4.1s) [68]. We used annually neutral net terrestrial exchange fluxes and turned off ship and aviation emissions and the 3-D Chemical Oxidized Source. Fossil fuel CO₂ emissions used in these tracer models were from the Open-source Data Inventory for Anthropogenic CO₂ (ODIAC), which is a global anthropogenic CO₂ emission dataset based on estimates made by the Carbon Dioxide Information and Analysis Center (CDIAC) and satellite-observed nightlights and the power plant profile dataset for spatial and temporal distributions [51]. It is monthly data with a spatial resolution of 1 × 1 degree.

The global background and a nested GEOS-Chem simulated XCO₂ for Eastern Asia (referred as Back-GEOS-XCO₂ and Nested-GEOS-XCO₂, respectively) were used for local XCO₂ background corrections without introducing modeled local biosphere flux variability. In the nested Eastern Asia CO₂ simulation, we used the GFED 4.1s land net ecosystem exchange, CarbonTracker biosphere prior CO₂ fluxes with annually neutral net CO₂ fluxes instead of the CarbonTracker biosphere posterior CO₂ fluxes, which are time-varying from year-to-year. The boundary conditions used for the GEOS-Chem nested grid model in eastern Asia at 0.5° (latitude) × 0.625° (longitude) simulation were provided by the global background model at 2° (latitude) × 2.5° (longitude) horizontal resolution. We made a restart file with 3-D CO₂ concentrations from CarbonTracker at 0:00, on Jan 1 of 2003. The resolution from 25 vertical layers and 2° × 3° in CT to 47 vertical layers and 2° × 2.5° in GEOS-Chem, was resampled using the nearest neighbor method. Model CO₂ profiles (averages for 12:00–14:00 local time) and the pressure-weighting function described in [66] were used to convert modelled CO₂ in discrete layers to XCO₂. Temporal mean values of GEOS-Chem background and GEOS-Chem nested simulated XCO₂ from 2003 to 2016 over Eastern Asia are shown in Figure S6.

In addition, changes in the sum of XCO₂ components (XCO₂ changes resulting from fossil, fire, biosphere, and ocean CO₂ fluxes) from CT and the MACC-III inversion models were used to estimate relationships between local surface CO₂ fluxes and XCO₂ variability described in Section 2.3.

2.2.4. XCO₂ Detrending Methods

To quantify local biosphere-driven changes in XCO₂ over time, we need to remove the background XCO₂ increase due to sustained fossil fuel emissions and influence from ocean exchange and fire, as well as remote biospheric influence from the study region [69,70]. The precise influence of all non-local non-biospheric processes cannot be directly measured, so we tested four different approaches for defining the background trend to isolate the local biospheric influence in the GM-XCO₂ record. For a first approximation, we adopted the global mean yearly increase from NOAA (https://www.esrl.noaa.gov/gmd/ccgg/trends/gr.html) to build a global background XCO₂ increase and applied that uniformly to all grid cells in the study region. The three other detrending methods were based on curve-fitting long-term trends, calculated and applied on a grid cell-by-cell basis. For one, we used the long-term trend XCO₂ increase from the gridded GM-XCO₂ observations themselves. In another, we used the fitted XCO₂ gridded results from CT-XCO₂. Finally, we fitted the global background and nested GEOS-Chem simulated XCO₂ to clarify the remote influences from the local influences.

The curve-fitting equation applied to each of the three XCO₂ records used a combination of a step-wise linear function and annual periodic function as shown in Equation (4).

$$f\left(t\right)=k_0+k_{1}t+{{\displaystyle \sum }}_{m=1}^4\left(a_i\cos \left(2\pi mt\right)+b_i\sin \left(2\pi mt\right)\right)$$

(4)

Here, $f\left(t\right)$ is XCO₂ in time t, $k_0$ and $k_{1}t$ are treated as the background and cumulated increase of XCO₂ in each year, and ${{\displaystyle \sum }}_{i=1}^4\left(a_i\cos \left(2\pi mt\right)+b_i\sin \left(2\pi mt\right)\right)$ is used for characterizing the XCO₂ seasonal cycle. An example of GM-XCO₂ and CT-XCO₂ fitting of 13 years in one grid cell (Lat: 38.5°N; Lon: 115.5°E) over EC is shown in Figure 3. Because fitting trends for individual years leads to discontinuities due to the start and end values, we used the mean value of three fitting results using 3 different time windows, starting from 1 December of previous year, 1 January, or 1 February of that year. Example curve fits are shown in Figure 3, but different patterns may be observed at other grid cells. The background XCO₂ increase (y) in year i and day j are extracted from the non-seasonal portion of Equation (4) and can be interpolated to daily values using Equation (5).

$$y_{i,j}=k_{i0}+{\mathrm{k}}_{\mathrm{i}1}\ast \mathrm{j}/n$$

(5)

where $k_{i0}$ and $k_{i1}$ are the fitted XCO₂ background ($k_0$) and increase ratio ($k_1$) in year i (from Equation (4) fit), n is the total units in one year (‘365′ used for daily data and ‘46′ used for data with time resolution of 8 days), and j is the current time step (day or number of 8-day interval). The yearly fitting performed well in global/site CO₂ yearly increase with differences of −0.04 ± 0.27 ppm for NOAA global mean, −0.03 ± 0.13 ppm for global mean of CT-XCO₂, and −0.06 ± 0.19 ppm for Mauna Loa shown in Figure S7 and Table S1.

2.2.5. Quantifying XCO₂ Trends

Global detrended GM-XCO₂ maps derived from integrated and gap-filled satellite observations with temporal and spatial resolution of 8 days and 1 × 1 degree were examined for growing period (GP) drawdown trends from 2003 to 2015. The GP was defined for each grid cell using the EVI analysis presented in Section 2.1.2 and the mean GP value of the detrended GM-XCO₂ was calculated for each grid cell and for each year. Least-squares linear fitting and the Mann–Kendal test were used to calculate XCO₂ trends and significance for each grid cell.

2.3. Estimating Surface CO₂ Fluxes from XCO₂ Variability Based on Correlations in Inversion Models

XCO₂ variability is closely related to surface CO₂ fluxes, however quantifying that relationship is challenging due to uncertainties in atmospheric transport and remote versus local CO₂ contributions. To estimate the magnitude of the surface flux responsible for observed trends in XCO₂, we derived a relationship between surface CO₂ fluxes and atmospheric XCO₂ from model simulations including CarbonTracker and MACC-III across different time and spatial scales: Daily and 16 days and global and grid cells in EC regions (shown in Figure 1). The total area of selected grid cells in EC shown in Figure S8 is 4.47 × 10¹² m². As the regions become smaller, we expected more noise in the correlations due to transport of CO₂ in/out of the region. We also examined the relationship over different time intervals. Daily XCO₂ changes were calculated from the difference between two consecutive days and compared to the daily mean CO₂ flux. The daily XCO₂ and CO₂ flux changes from 1 January 2003 to 31 December 2016 in each grid cell were calculated using CT and MACC. XCO₂ changes of several days were defined as the sum of daily changes and CO₂ flux changes were defined as total CO₂ fluxes of these days. Correlation analysis and least-squares fitting were calculated for the flux-XCO₂ relationships (defined as transform ratios) at these different time and space scales. For this analysis of EC, we chose the 16-day correlation over EC to estimate CO₂ flux anomalies from XCO₂ changes.

2.4. Supporting Independent Datatsets

2.4.1. Vegetation Coverage Changes

We compared trends in atmospheric XCO₂ to independent evidence of vegetation changes in EC. We adopted tree canopy cover percentage (TCCP) from global vegetation continuous fields (VCF) products by Song et al. [50], which were retrieved from the NASA website (https://search.earthdata.nasa.gov/). This vegetation continuous field was developed using daily satellite AVHRR observations. Tree canopy cover percentage with a spatial resolution of 0.05 × 0.05 degree was resampled to 1 × 1 degree, same as that of vegetation parameters and GM-XCO₂ used in this analysis.

2.4.2. XCO₂ from TCCON Measurements

TCCON is a global network of ground-based Fourier transform spectrometers for atmospheric XCO₂ and other trace gases measurements [71]. We compared GM-XCO₂ retrievals to these ground-based measurements to identify potential measurement biases. XCO₂ retrievals from TCCON has a high accuracy with approximately 0.25%, which can be used for validation of satellites observed XCO₂ [72]. In this study, we selected atmospheric CO₂ measurement from TCCON sites with data between 2003.01 to 2016.03 for GM-XCO₂ validation. These sites include Bialystok, Bremen, Karlsruhe, Orleans, Garmisch, Park Falls, Lamont, Tsukuba, Edwards, JPL/Caltech, Saga, Manus, Darwin, Wollongong, and Lauder. Locations of former 11 sites in northern hemisphere are shown in Figure 1, and other the 4 sites (Manus, Darwin, Wollongong, and Lauder) in the southern hemisphere are inside of Brazil, north and south of Australia, and New Zealand, respectively. TCCON measurements from Jan 2003 to Mar 2016, observing local time during 11:00 to 15:00 were used for XCO₂ comparison. As a result, there were 3603 pairs of GM-XCO₂ and TCCON XCO₂, from 15 TCCON sites, used for comparison (GM-XCO₂ minus TCCON XCO₂). R², mean bias, mean absolute bias, and standard deviation between them were 0.96, 0.20, 0.90, and 1.16, respectively.

2.4.3. Climate Factors

Land surface temperature (LST) and precipitation were used for evaluating the effects from climate change in XCO₂ change. LST from MODIS (MOD11C2 v6) [52] was downloaded from NASA (https://search.earthdata.nasa.gov/). Precipitation from the Tropical Rainfall Measuring Mission (TRMM) v3B42 [53] was obtained from Google Earth Engine. LST and precipitation were resampled to temporal and spatial resolution of 16 days and 1 × 1 degree for correlation analysis with corrected XCO₂.

3. Results

3.1. Trends in Regional Vegetation Parameters

Here, we show mean growing period (GP) vegetation trends for MODIS NDVI and AVHRR NDVI from 2003 to 2016 and AVHRR LAI and tree canopy cover trends from 2003 to 2016 (Figure 4). These vegetation indices are unitless because NDVI and EVI range 0–1 and LAI is in m² of leaf area per m² land surface area. MODIS EVI, LAI, and GPP trends during the GP from 2003 to 2016 are shown in Figure S9. Significant greening trends were observed in EC (Eastern China with red edge) with similar patterns in MODIS NDVI, AVHRR NDVI, and MODIS EVI (Figure 4a,b and Figure S9a). The trends in mean GP values of MODIS NDVI, AVHRR NDVI, and MODIS EVI over EC from 2003 to 2016 were 0.38 ± 0.27 × 10⁻² yr⁻¹, 0.38 ± 0.30 × 10⁻² yr⁻¹, and 0.32 ± 0.21 × 10⁻² yr⁻¹, which accounted for 0.70 ± 0.61% yr⁻¹ and 0.71 ± 0.75% yr⁻¹ of NDVI (MODIS and AVHRR, respectively) and 0.86 ± 0.60% yr⁻¹ of EVI relative to the mean values over this 14-year period. The standard deviation reported here is the variability in all grid cells trends across the EC domain and not an indicator of within-grid cell trend uncertainty. AVHRR LAI also showed significant increasing trends in EC with 0.013 ± 0.010 m² m⁻² yr⁻¹ (Figure 4c). In addition, increasing trends (0.009 ± 0.010 m² m⁻² yr⁻¹ and 0.11 ± 0.19 g C m⁻² yr⁻²) of LAI and GPP from MODIS were also significant in most area of Eastern China (Figure S9). Total increases of NDVI and GPP over Eastern China from 2003 to 2015 could be around 8.40% to 8.52% and 9.24 g C m⁻², respectively.

Tree canopy annual trends from 2003 to 2016 showed a strong net increasing trend (0.65 ± 0.34% yr⁻¹) in the southern part of EC (Figure 4d, latitude < 34°N). Large total tree canopy increase from 2003 to 2015 (7.80%), potentially resulting in enhanced photosynthesis rates, may have increased the net carbon sink in this region.

3.2. Trends in Regional XCO₂ Summer Drawdown

Here we identify the trends in GM-XCO₂ that best reflect local biospheric processes in EC and present different methods for removing the influence from fossil fuel and fire emissions and ocean exchange from the observations, and in some cases, remote biospheric influence.

3.2.1. Satellite Observed XCO₂ Trends after Removing Global Mean Growth

As a first attempt, we removed the NOAA global mean XCO₂ growth rate from the GM-XCO₂ in order to calculate summer drawdown trends using the mean value of detrended GM-XCO₂ during the GP of each year as shown in Figure 5a. This removed remote influence, the long-term atmospheric growth rate, at each grid cell. Red grid cells with black outlines, indicating significant decreasing trends, were observed in some parts of south China (−0.07 ± 0.05 ppm yr⁻¹ for decreasing grid cells). Blue grid cells with black outlines, indicating significant increasing trends, 0.05 ± 0.05 ppm yr⁻¹, were observed in some urban regions like Beijing and Shanghai.

Regions around cites like Beijing and Shanghai showed increasing XCO₂ at more extreme rates than the global mean. Blue grid cells indicate CO₂ emission increases (0.045 ± 0.055 kg C m⁻² yr⁻²) near urban centers of Eastern Asia in Figure 5a. This is because fossil emissions have increased faster in Beijing and Shanghai than most other regions. Gridded trends in anthropogenic emissions from January 2003 to December 2016 from ODIAC are shown in Figure 5b. Therefore, it is necessary to remove the influence of the local urban CO₂ emission domes in order to isolate trends in local biospheric influence on XCO₂.

3.2.2. Satellite-Observed XCO₂ Trends after Removing Influence of the Urban CO₂ Dome

In order to remove the influence of increasing anthropogenic emissions from urban centers in our study area, we tested four different spatially variable methods of detrending GM-XCO₂ observations: GM-XCO₂ fitting (Figure 6a), CT-XCO₂ fitting (Figure 6b), and GEOS-Chem simulated XCO₂ fitting at two different spatial scales: a coarse 2 × 2.5 degree lat/lon simulation and a high-resolution nested simulation for Eastern Asia (Figure 6c–d). All methods showed GM-XCO₂ residuals decreasing in many EC grid cells from 2003 to 2015 indicating greater CO₂ uptake during the GP. For the entire EC domain, detrending the XCO₂ time series of each grid by the long-term trend of each grid cell (GM-XCO₂ detrending using the full record including growing season and dormant season) resulted in a small trend in the corrected GM-XCO₂ during the GP, −0.058 ± 0.029 ppm yr⁻¹. The standard deviation reported here is the variability in individual grid cell trends across the EC domain and not an indicator of within grid cell trend uncertainty. This is possible because the annual means and the GP drawdown only can have different trends as seasonal cycles change over time. Detrending using model simulations of XCO₂ from remote and local non-biospheric CO₂ fluxes resulted in larger GP CO₂ uptake trends: −0.084 ± 0.090, −0.081 ± 0.088, and −0.070 ± 0.093 ppm yr⁻¹ (mean ± 1 standard deviation in space), for CT-XCO₂ and GEOS-Chem coarse and fine resolution, respectively. Of all the grid cells, larger GP CO₂ uptake trends were observed in 80%, 79%, and 77% of the grid cells depending on the correction method applied, CT-XCO₂ fitting, and GEOS-Chem fittings at coarse and fine resolution, respectively. Grid cells with statistically significant (p < 0.05) increasing GP CO₂ uptake trends occupied 40%, 41%, and 39% of EC, respectively. The trends for just the grid cells that showed increased CO₂ uptake were −0.116 ± 0.067, −0.113 ± 0.069, and −0.108 ± 0.068 ppm yr⁻¹, respectively (mean ± 1 standard deviation in space).

The similar pattern of trends of these last three correction methods showed that there was robust evidence of increased GP CO₂ uptake in EC due to biological activity. We consider the corrected GP trend using the GM-XCO₂ annual increase of −0.058 ppm yr⁻¹ to be the lower estimate of biospheric influence on EC CO₂ uptake because the local biospheric influence is in the GM-XCO₂ signal itself and a portion of it is removed by the fitting. CT-XCO₂ should depict the general pattern of atmospheric CO₂ concentration over city areas, as the posterior CO₂ concentration and flux from CT includes fossil fuel emissions and was the optimum estimate of the inversion model.

The spatial patterns of corrected GM-XCO₂ decreases were similar to the increases in vegetation indices, south of the urban centers, shown in Figure 4 and Figure S9a. These results indicate vegetation increases have been accompanied by stronger photosynthesis rates during the GP from 2003 to 2016 and imply increased carbon uptake in EC during the growing season.

3.3. Trends in Flux Estimates

The corrected GM-XCO₂ trends of −0.070 to −0.084 ppm yr⁻¹ in EC suggest an increased biospheric uptake. In order to estimate the magnitude of the surface CO₂ fluxes that would produce that XCO₂ anomaly, we used the CT and MACC models to relate CO₂ fluxes to atmospheric CO₂ concentration anomalies. We explored correlation between model simulated XCO₂ changes daily or over 16 days and net CO₂ flux changes in the global or EC region shown in Figure 7. Significant and strong positive correlations were found in the global daily mean data both from CT (Figure 7a, R² = 0.84, p < 0.01) and MACC (Figure 7d, R² = 0.49, p < 0.01) with slopes approaching 0.26 ppm/(g C/m²/day) from CT and 0.21 ppm/(g C/m²/day) from MACC, where the former was a little larger. The correlation and slopes over EC were much lower than that in the global mean, but still significant (Figure 7b–d). This is because at the global scale, CO₂ is conserved during atmospheric transport within the global domain, whereas at the regional domain, transport can be a source or sink of CO₂. The correlation decreased as the domain considered got smaller from global to EC. The correlation from daily to 16 days (the time resolution of satellite observed XCO₂ and vegetation parameters) became stronger both for CT (Figure 7c) (R² = 0.26, p < 0.01) and MACC (Figure 7f) (R² = 0.23, p < 0.01). The flux-XCO₂ transform ratios of 16-day were similar for CT (0.0523 ppm/(g C/m²/16days)) and MACC (0.0465 ppm/(g C/m²/16days)). The correlation became stronger as the time interval increased (from daily to 16-day means), because the influence of persistent CO₂ fluxes built up within the domain.

We estimated the CO₂ flux change from the corrected GM-XCO₂ trends and the model correlation between XCO₂ 16-day changes and CO₂ flux change shown in

Table 2. GM-XCO₂ decreasing trends and estimated CO₂ uptake changes from 2003 to 2015 during the vegetation growing period over Eastern China (EC). There are results from three different XCO₂ background corrections and two flux-XCO₂ transform ratios. XCO₂ background corrections are shown in Section 3.2. Transform ratios are shown in Figure 7. Uncertainties in XCO₂ trends and changes refer to 1 standard deviation of spatial grid-scale variability in EC.

Attributes\Corrections	CT-XCO2		Back-GEOS-XCO2		Nested-GEOS-XCO2
XCO2 trends (ppm/yr)	−0.084 ± 0.090		−0.081 ± 0.088		−0.07 ± 0.093
XCO2 change (ppm)	−1.008 ± 1.080		−0.972 ± 1.056		−0.840 ± 1.116
Transform ratios	CT	MACC	CT	MACC	CT	MACC
Carbon uptake increase (g C/m2)	9.57 ± 10.94	10.30 ±11.85	8.88 ± 10.48	9.52 ± 11.33	6.36 ± 11.63	6.69 ± 12.62
Carbon amount (Tg C)	42.76	46.04	39.69	42.58	28.41	29.89

. The total XCO₂ change were estimated as −1.008 ± 1.080, −0.972 ± 1.056, and −0.840 ± 1.116 ppm (standard deviation represents spatial variability) from 2003 to 2015 for different correction methods (CT-XCO₂ fitting, background GEOS-XCO₂ fitting, and nested GEOS-XCO₂ fitting). The inferred carbon uptake increase mean value over the region and its spatial variability were 9.57 ± 10.94, 8.88 ± 10.48, and 6.36 ± 11.63 g C/m² using the CT ppm-to-flux transform ratio and 10.30 ±11.85, 9.52 ± 11.33, and 6.69 ± 12.62 g C/m² using the MACC transform ratio. That is basically consistent with GPP increase of 9.24 g C m⁻² shown in Section 3.1. The total carbon uptake amount from 2003 to 2015 ranged from 28.41 to 46.04 Tg C with an EC region area of 4.47 × 10¹² m² shown in Figure S8.

To put this carbon uptake into perspective, we compared it to the carbon exchange inferred from CT during the GP (from day 126 to 275) and the whole year in this region. The averaged biospheric CO₂ flux during the GP during 2003 to 2016 over EC was −99.8 ± 30.6 g C m⁻² yr⁻¹ (uncertainty is 1 standard deviation of inter-annual variability). Annual fossil fuel CO₂ emission for the same period was 328.48 ± 91.38 g C m⁻² yr⁻¹ (uncertainty is 1 standard deviation of inter-annual variability). The satellite-observed total carbon uptake increases from 2003 to 2015 accounted for 6.4% to 10.3% of the mean biosphere flux during the GP, approaching the percentage increase in the vegetation indices of 8.4% to 10.3% from 2003 to 2015. This uptake could offset 1.9% to 3.1% of annual fossil carbon emissions over that region.

3.4. Possible Contribution from Climate Variability, Independent of Land Management

To check the relationship between XCO₂ and climate change, we present the Pearson correlation coefficients between corrected XCO₂ and MODIS land surface temperature (LST) and TRMM precipitation during the GP from 2003 to 2016 shown in Figure 8. In EC, there was no significant correlation between corrected XCO₂ and temperature. This is because vegetation growth is not temperature-limited due to already warm temperatures (Figure S10a) and there were no significant grid-cell trends in LST (Figure S10c). The correlation between corrected XCO₂ and precipitation was also not significant. The precipitation in the southern part of EC increased and in the northern part decreased (Figure S10d). A decrease in rainfall in that region was unlikely to result in increased CO₂ uptake. Vegetation growth was limited by water in the northern region (Figure S10b), which is why water diversion projects have been constructed to transfer water to the north, therefore disconnecting CO₂ uptake from precipitation. Overall, there is no strong support for the observed EC greening trends to be caused by changes in regional temperature and rainfall. It is more likely that greening patterns are driven by increased irrigation in Northeastern China [22,73] and land use change in Southeastern China [74].

4. Discussion

Vegetation greening trends (Figure 4) and XCO₂ decreasing trends (Figure 6) during the GP show similar patterns of enhanced greening and CO₂ uptake in Eastern China (EC). This indicates that enhancements in vegetation activity have had a measurable influence on regional atmospheric XCO₂ concentrations due to enhanced CO₂ uptake during the GP. Here, we discuss the likely drivers of the vegetation changes we observed, our estimates of the flux changes associated with the XCO₂ trend, and consider the uncertainty in our analysis.

4.1. Different Drivers of Greening in North and South Parts of Eastern China

Greening trends are observed throughout EC (Figure 4), however the main drivers of these trends are likely different in the north and the south. In the northern part, cropland production has increased significantly from 1987 to 2010 [10] in order to feed the large increase in Chinese population (from 1.08 to 1.34 billion). Increases in agricultural productivity come from increased cropland area [10,74,75] and increased irrigation from water diversion projects moving water from the south to the North China Plain [22,73]. In addition, increased fertilization by nitrogen and phosphorus application [76] and lengthening of the growing period due to warmer temperatures with climate change [77] could provide more nutrients and a longer period of optimal growth conditions.

In the southern part, tree canopy coverage has increased by 9.27 ± 4.93% in 2016 from 30.51 ± 13.90% in 2003 [50]. Over the last century, widespread deforestation and forest degradation occurred due to conversion to cropland, infrastructure development, mining, fire, and urbanization. During the 21st century, afforestation and reforestation increased significantly in the south part of EC [4], due to governmental policy. Climate conditions over southern EC are suitable for forest growth and expansion to increase carbon uptake with protection provided by land management policies [9].

4.2. Uncertainty of GM-XCO₂ Trends

The GM-XCO₂ product involves merging and gap-filling multiple satellite records (Section 2.2.2), therefore, we need to consider the uncertainty associated with this product. Spurious trends could be caused by the gap-filling method, so we checked the amount of available satellite-observed XCO₂ in each grid cell over EC (Figure S11). A greater proportion of observations are distributed in the northern part of China for SCIAMACHY BESD XCO₂ (2003.01–2009.05), GOSAT ACOS XCO₂ (2009.06–2014.08), and OCO-2 XCO₂ (2014.09–2016.03) because of the cloud contamination discussed earlier. Here, we show the temporal change of integrated XCO₂ (the merging of data from the three sensors) and GM-XCO₂ (gap-filled) using available data in EC as shown in Figure 9. GM-XCO₂ agreed well with integrated XCO₂ with a mean difference of 0.25 ± 1.33 ppm. Some differences exist during summer for distinct XCO₂ spatial variations shown with mean ± 1 STD GM-XCO₂ (red dashed lines). We found no significant trend in the GP residuals between integrated XCO₂ and the gap-filled GM-XCO₂ using the Mann–Kendal trend test. Therefore, we conclude that the XCO₂ decreasing trends are not caused by the gap-filling method.

In addition, we checked the GM-XCO₂ product against TCCON measurements as shown in Figure S5. The mean absolute difference (bias) between GM-XCO₂ and TCCON observations with 3608 co-located data pairs was 0.90 ± 1.16 ppm. Individual grid-cell GM-XCO₂ bias followed a Gaussian distribution, therefore the XCO₂ uncertainty over EC was estimated as the averaged bias divided by the square root of the number of grid cells in EC. As a result, the uncertainty was less than 0.045 ppm, which is significantly less than the XCO₂ total change (0.840–1.008 ppm from 2003 to 2016). XCO₂ differences between GM-XCO₂ and TCCON measurements from 2003 to 2016 are shown in Figure 10. The differences fluctuate around zero and present no trends from 2003 to 2016 for data pairs in the northern hemisphere. This comparison with the TCCON observations showed that there was no evidence of sensor offsets or drift that would create artificial GM-XCO₂ trends over time.

4.3. Uncertainties in Constraining the Urban Dome and Remote Influence on XCO₂ Trends

The XCO₂ trends observed over EC result from a combination of remote and local influences. Accurately quantifying the influence of regional biospheric CO₂ fluxes in EC depends on how well we remove the influence from non-local biospheric CO₂ fluxes including remote and local anthropogenic fossil and fire emissions and ocean gas exchange as well as remote biospheric influence. In this study, we tested four methods for removing the cumulated XCO₂ increase to better study the local biospheric component of XCO₂ changes. The first method using the NOAA-reported global mean CO₂ is not able to remove the influence from the local urban CO₂ dome, which increases over time and therefore underestimated the trend in XCO₂ due to local biospheric influence. In addition, the time series fitting of GM-XCO₂ removed the local annual biospheric contribution to the XCO₂ trend and underestimated the signal of interest, although there was still a component of GS trends in XCO₂ that served as a lower bound of local biospheric influence. The trends calculated from detrending using simulated XCO₂ using CT and GEOS-Chem transport (with and without nested grid cells) of CT surface CO₂ fluxes were likely the best representation of local biospheric influence on XCO₂ because they are designed to remove trends caused by the remote biospheric CO₂ fluxes over time, but retain trends in local biospheric CO₂ fluxes, plus remove all the non-biospheric influence. Similarities between Figure 6c,d showed that characterizing the fine spatial scale resolution of the urban XCO₂ dome was not critical to defining GP trends. Therefore, XCO₂ decreasing trends of −0.084 ± 0.090, −0.081 ± 0.088, and −0.070 ± 0.093 ppm yr⁻¹ were our best estimates of XCO₂ change caused by EC region local biospheric CO₂ fluxes with uncertainty listed is the grid-scale spatial variability in trends (Table 2). These trends result in a total XCO₂ change from 2003 to 2015 of −0.84 to −1.01 ppm.

This analysis is subject to two sources of bias. First, uncertainty in remote CO₂ fluxes could result in artificial trends in local biospheric CO₂ fluxes in the study region, however it is very unlikely that such an artificial signal from remote CO₂ flux uncertainty would correspond to the exact region that shows vegetation greening. Second, an overestimate in local fossil fuel CO₂ emissions from EC in the EDGAR database would result in an overestimate of local biospheric uptake. Saeki and Patra [78] suggested that fossil fuel CO₂ emissions in EC may in fact be overestimated based on scaling with methane emissions. There is no independent way to verify local fossil emissions in EC at the moment. However, the correction method based on detrending using annual global mean CO₂ concentrations from NOAA supports a local biospheric sink in the southern portion of EC in addition to showing urban emissions increasing faster than the global mean in the northern part (Figure 4a). This correction method cannot create regional artificial biospheric uptake due to local CO₂ emissions uncertainty. However, the magnitude of the inferred sink is sensitive to local CO₂ emissions uncertainty. This demonstrates the importance of independent constraints of fossil and biospheric CO₂ fluxes in order to understand CO₂ flux trends near urban centers.

4.4. Inferred Surface Flux Changes over This Period

It is predicted that China’s investments into environmental improvement has likely increased the forest carbon sink in EC [9,23,24]. However, the carbon sink over the East Asian monsoon region is difficult to measure and there is speculation that it has been underestimated [79]. We find the estimated cumulated carbon uptake increases from 2003 to 2015 was in range from 6.36 to 10.30 g C m⁻² over EC. To understand the significance of these changes, we compare with mean posterior flux estimates from CT (

Table 3. Comparison of the estimated carbon uptake increase (6.36 to 10.30 g C m⁻²) due to greening with various posterior CO₂ fluxes inferred from CT 2017 for the region during the growing period (GP) and the whole year. Posterior CO₂ fluxes include the net flux, biosphere (Bio) flux, and fossil fuel emissions of annual mean and inter-annual variation. Negative values indicate carbon uptake and positive values indicate carbon emission. Values in parentheses compare the CO₂ sink enhancement as percentages of the various GP CT CO₂ fluxes.

Period	Net Flux (g C m−2 yr−1)	Bio Flux (g C m−2 yr−1)	Fossil Flux (g C m−2 yr−1)
GP	28.0 ± 54.5	−99.8 ± 30.6	125.7 ± 34.8
	(22.6% to 36.7%)	(6.4% to 10.3%)	(5.1% to 8.2%)
Annual	264.15 ± 109.85	−70.77 ± 36.76	328.48 ± 91.38
	(2.4% to 3.9%)	(9.0% to 14.6%)	(1.9% to 3.1%)

). The carbon sink increase is equivalent to 9.0% to 14.6% of the annual biospheric flux over EC and equivalent to offsetting about 1.9–3.1% of EC’s annual fossil fuel CO₂ emissions. By comparison, the increased carbon uptake during the GP was slightly underestimated in the interannually varying CT posterior CO₂ fluxes, which is about 6.4% to 10.3% of regional CT posterior annual biopsheric flux.

In this estimation, we calculated the biospheric CO₂ flux trend during the GP only, not including the dormant season. It is possible that carbon release during the dormant season could cancel some of the carbon uptake gains made in the GP. We tested XCO₂ trends in the dormant season, after CT-XCO₂ correction (Figure S12), but find no significant trends in EC. It is possible that this carbon sink could saturate in the future as the forests mature and agricultural intensification peaks [80]. These could be evaluated with continuous observations.

The method used here to infer flux changes is different from atmospheric inversion approaches [26,27] and more recent Lagrangian methods [81,82] to link surface CO₂ fluxes to atmospheric CO₂ observations. Our approach may likely underestimate the surface CO₂ sink strength as some of the negative CO₂ anomaly is transported out of the study domain. However, it is unlikely that an anomaly caused by a CO₂ sink outside of the domain is transported inside because of the distribution of vegetation in this region. This analysis should be viewed as an independent method with a different sensitivity to atmospheric transport uncertainty.

5. Conclusions

Significant greening trends in Eastern China were observed by various parameters including NDVI trends of 0.70 ± 0.61% yr⁻¹ (MODIS) and 0.71 ± 0.75% yr⁻¹ (AVHRR) and an EVI trend of 0.86 ± 0.60% yr⁻¹ (MODIS). Greening trends over the northern and southern parts of EC are consistent with increased cropland area and production increase in north, and forest canopy increase in south. These have previously been attributed to the national policies of sustainable development, including intensification of agriculture and afforestation and reforestation as part of China’s Grain-to-Green program [2]. Here, we show that these greening patterns are consistent with increased growing season CO₂ uptake in this region.

Satellite-observed XCO₂ trends showed a corresponding increase in the carbon sink in the region of greening in EC. Significant XCO₂ decreases attributable to local biospheric influence during the growing period were shown using different correction methods to remove non-local, non-biospheric influence. Averaged XCO₂ trends across EC varied from −0.070 to −0.084 ppm yr⁻¹. Estimated local CO₂ uptake increases were 6.36 to 10.30 g C m⁻² over EC during the growing period, which is consistent with a MODIS observed GPP increase of 9.24 g C m⁻². These fluxes accounted for 6.4 to 10.3% of the mean biosphere flux during the GP, which is similar to percent increases in vegetation indices, 8.40 to 10.32% from 2003 to 2015. The total carbon uptake increase ranged from 28.41 to 46.04 Tg C. That could offset about 2–3% of EC’s annual fossil fuel emissions. Integrated data from multiple greenhouse satellite (SCIAMACHY, GOSAT, OCO-2 or the recently launched Tansat) observations could provide a valuable method for estimating regional CO₂ flux change.