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

: Persistent and widespread increase of vegetation cover, identiﬁed 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 di ﬃ cult 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 2 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 2 dry air mole fraction (XCO 2 ) from January 2003 to May 2016, based on observations from SCIAMACHY, GOSAT, and OCO-2 satellites, were used to calculate XCO 2 changes during growing season for 13 years. We derived a relationship between XCO 2 and surface net CO 2 fluxes from two inversion model simulations, CarbonTracker and Monitoring Atmospheric Composition and Climate (MACC), and used those relationships to estimate the biospheric CO 2 flux enhancement based on satellite observed XCO 2 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 2 emissions, we estimated an enhanced XCO 2 drawdown during the GP of − 0.070 to − 0.084 ppm yr − 1 . Increased carbon uptake during the GP was estimated to be 28.41 to 46.04 Tg C, mainly from land management, which could o ﬀ set about 2–3% of EC’s annual fossil fuel emissions. These results show the potential of using multi-satellite observed XCO 2 to estimate carbon ﬂuxes 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.


Introduction
Greening trends appear over much of the global land area [1], indicating changes in ecosystem function at regional scales. Atmospheric carbon dioxide (CO 2 ) 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 Datasets used in this research are summarized in Table 1. 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) XCO2 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) CO2 fluxes and corresponding modeled XCO2 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.  Datasets used in this research are summarized in Table 1. 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 2 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 2 fluxes and corresponding modeled XCO 2 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

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).
where t is time in days, y(t) 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).

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.

Satellite-Observed XCO2
Multiple satellite-observed XCO2 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 XCO2 was retrieved by the full-physics based Bremen Optimal Estimation-DOAS (BESD) algorithm (v02.01.01), referred to as BESD-XCO2 [43,44] and was downloaded from the European Space Agency (ESA, http://www.esa-ghg-cci.org/). GOSAT was retrieved by the Atmospheric CO2 Observations from Space (ACOS) project team (v7.3), referred to as ACOS-XCO2 [45]. OCO-2 was retrieved by the OCO-2 team (v7r), referred as OCO-XCO2 [46]. ACOS-XCO2 and OCO-XCO2 were downloaded from NASA (https://disc.gsfc.nasa.gov/). XCO2 from these different observations during 2003 to 2016 was integrated and mapped, described in Section 2.2.2, into global maps of XCO2 (referred as GM-XCO2) with temporal and spatial resolution of 8 days and 1 × 1 degree. for all years, mean seasonal EVI (red points), and EVI fitting results (black for spring increase and blue for fall decrease) are shown in the top panel. Curvature and derivative of the curvature of EVI fitting (black for increase and blue for decrease) are shown in bottom panel. The period between the two red dash lines was identified as the growing period for this grid cell.

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. Multiple satellite-observed XCO 2 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 2 was retrieved by the full-physics based Bremen Optimal Estimation-DOAS (BESD) algorithm (v02.01.01), referred to as BESD-XCO 2 [43,44] and was downloaded from the European Space Agency (ESA, http://www.esa-ghg-cci.org/). GOSAT was retrieved by the Atmospheric CO 2 Observations from Space (ACOS) project team (v7.3), referred to as ACOS-XCO 2 [45]. OCO-2 was retrieved by the OCO-2 team (v7r), referred as OCO-XCO 2 [46]. ACOS-XCO 2 and OCO-XCO 2 were downloaded from NASA (https://disc.gsfc.nasa.gov/). XCO 2 from these different observations during 2003 to 2016 was integrated and mapped, described in Section 2.2.2, into global maps of XCO 2 (referred as GM-XCO 2 ) with temporal and spatial resolution of 8 days and 1 × 1 degree.

Merging and Gap-Filling XCO 2 Data Products (Integration and Mapping)
In order to construct a long-term XCO 2 dataset, BESD-XCO 2 from January 2003 to May 2009,  ACOS-XCO 2 from June 2009 to August 2014, and OCO-XCO 2 from Sep 2014 to March 2016 were used for XCO 2 integration and mapping. XCO 2 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 2 data, and (3) global land mapping of XCO 2 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 2 from CarbonTracker (CT-XCO 2 ) selected as the common proxy for adjustment (Equation (2)) [43,62].
where XCO 2adj,t is the adjusted XCO 2 value at observation time t, XCO 2ret,t corresponds to retrieved XCO 2 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 2 profile from CT, and X a,t is a priori CO 2 profile from satellite retrievals.
Considering the diurnal variability in XCO 2 and differences in satellite overpass times, step 2 unified XCO 2 from multiple satellites to a common observing local time and spatial resolution. We introduced a scale factor derived from the diurnal variation of XCO 2 using CT-XCO 2 to normalize the overpass time of XCO 2 retrievals from these three satellites to a reference time (Equation (3)).
where XCO 2nor,t is the normalized satellite-observed XCO 2 at the reference time (rt: 13:00, local time); XCO 2adj,t is the adjusted satellite XCO 2 derived from Equation (2) at satellite overpass time t; XCO 2CT,rt and XCO 2CT,t are the XCO 2 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 XCO 2nor,t of BESD-XCO 2 , ACOS-XCO 2 , and OCO-XCO 2 respectively. This uniform spatio-temporal scale was introduced to aggregate different XCO 2 data products in order to maintain as much as possible the spatial (regional and national) and temporal (monthly and seasonally) variation of XCO 2 observed by these three satellites. The integrated XCO 2 dataset in 30 × 30 km and 8 days was generated from 2003 to 2016, hereafter referred to as integrated-XCO 2 . The data quantity and trends of integrated-XCO 2 from 2003 to 2016 is shown in Figure S2.
In step 3, we applied a gap-filling interpolation method to the integrated XCO 2 to produce continuous global XCO 2 maps in time and space. The spatio-temporal geostatistical mapping method, GM-XCO 2 , 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 2 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 2 from Europe and Asia was used in the XCO 2 trend analysis over EC.

Model Simulations of XCO 2 and CO 2 Fluxes
Estimated surface CO 2 fluxes and XCO 2 concentrations from the GEOS-Chem, CarbonTracker, and MACC models were used for two purposes: (1) To create estimates of XCO 2 background increase in the study area (Section 2.2.4), and (2) estimate relationships between local surface CO 2 fluxes and XCO 2 variability (Section 2.3). CarbonTracker (CT) 2017 is a global inversion model from NOAA providing atmospheric CO 2 mixing ratio distributions and inferred surface CO 2 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 2 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 2 total files. Contributions to XCO 2 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 2 and net CO 2 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 2 flux influence in Eastern China from remote biospheric and other non-biospheric CO 2 flux influences on that region, we created estimates of various CO 2 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 2 simulation, we used the CarbonTracker fossil fuel, ocean exchange, and biosphere posterior CO 2 fluxes to replace the standard GEOS-Chem surface CO 2 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 2 emissions used in these tracer models were from the Open-source Data Inventory for Anthropogenic CO 2 (ODIAC), which is a global anthropogenic CO 2 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 2 for Eastern Asia (referred as Back-GEOS-XCO 2 and Nested-GEOS-XCO 2 , respectively) were used for local XCO 2 background corrections without introducing modeled local biosphere flux variability. In the nested Eastern Asia CO 2 simulation, we used the GFED 4.1s land net ecosystem exchange, CarbonTracker biosphere prior CO 2 fluxes with annually neutral net CO 2 fluxes instead of the CarbonTracker biosphere posterior CO 2 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 2 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 2 profiles (averages for 12:00-14:00 local time) and the pressure-weighting function described in [66] were used to convert modelled CO 2 in discrete layers to XCO 2 . Temporal mean values of GEOS-Chem background and GEOS-Chem nested simulated XCO 2 from 2003 to 2016 over Eastern Asia are shown in Figure S6.
In addition, changes in the sum of XCO 2 components (XCO 2 changes resulting from fossil, fire, biosphere, and ocean CO 2 fluxes) from CT and the MACC-III inversion models were used to estimate relationships between local surface CO 2 fluxes and XCO 2 variability described in Section 2.3.

XCO 2 Detrending Methods
To quantify local biosphere-driven changes in XCO 2 over time, we need to remove the background XCO 2 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 2 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 2 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 2 increase from the gridded GM-XCO 2 observations themselves. In another, we used the fitted XCO 2 gridded results from CT-XCO 2 . Finally, we fitted the global background and nested GEOS-Chem simulated XCO 2 to clarify the remote influences from the local influences.
The curve-fitting equation applied to each of the three XCO 2 records used a combination of a step-wise linear function and annual periodic function as shown in Equation (4).
Here, f (t) is XCO 2 in time t, k 0 and k 1 t are treated as the background and cumulated increase of XCO 2 in each year, and 4 i=1 (a i cos(2πmt) + b i sin(2πmt)) is used for characterizing the XCO 2 seasonal cycle. An example of GM-XCO 2 and CT-XCO 2 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 2 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).
where k i0 and k i1 are the fitted XCO 2 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 2 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 2 , and −0.06 ± 0.19 ppm for Mauna Loa shown in Figure S7 and Table S1.

Estimating Surface CO2 Fluxes from XCO2 Variability Based on Correlations in Inversion Models
XCO2 variability is closely related to surface CO2 fluxes, however quantifying that relationship is challenging due to uncertainties in atmospheric transport and remote versus local CO2 contributions. To estimate the magnitude of the surface flux responsible for observed trends in XCO2, we derived a relationship between surface CO2 fluxes and atmospheric XCO2 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 12 m 2 . As the regions become smaller, we expected more noise in the correlations due to transport of CO2 in/out of the region. We also examined the relationship over different time intervals. Daily XCO2 changes were calculated from the difference between two consecutive days and compared to the daily mean CO2 flux. The daily XCO2 and CO2 flux changes from 1 January 2003 to 31 December 2016 in each grid cell were calculated using CT and MACC. XCO2 changes of several days were defined as the sum of daily changes and CO2 flux changes were defined as total CO2 fluxes of these days. Correlation analysis and least-squares fitting were calculated for the flux-XCO2 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 CO2 flux anomalies from XCO2 changes.

Vegetation Coverage Changes
We compared trends in atmospheric XCO2 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

Quantifying XCO 2 Trends
Global detrended GM-XCO 2 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 2 was calculated for each grid cell and for each year. Least-squares linear fitting and the Mann-Kendal test were used to calculate XCO 2 trends and significance for each grid cell.

Estimating Surface CO 2 Fluxes from XCO 2 Variability Based on Correlations in Inversion Models
XCO 2 variability is closely related to surface CO 2 fluxes, however quantifying that relationship is challenging due to uncertainties in atmospheric transport and remote versus local CO 2 contributions. To estimate the magnitude of the surface flux responsible for observed trends in XCO 2 , we derived a relationship between surface CO 2 fluxes and atmospheric XCO 2 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 12 m 2 . As the regions become smaller, we expected more noise in the correlations due to transport of CO 2 in/out of the region. We also examined the relationship over different time intervals. Daily XCO 2 changes were calculated from the difference between two consecutive days and compared to the daily mean CO 2 flux. The daily XCO 2 and CO 2 flux changes from 1 January 2003 to 31 December 2016 in each grid cell were calculated using CT and MACC. XCO 2 changes of several days were defined as the sum of daily changes and CO 2 flux changes were defined as total CO 2 fluxes of these days. Correlation analysis and least-squares fitting were calculated for the flux-XCO 2 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 2 flux anomalies from XCO 2 changes.

Vegetation Coverage Changes
We compared trends in atmospheric XCO 2 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 2 used in this analysis.

XCO 2 from TCCON Measurements
TCCON is a global network of ground-based Fourier transform spectrometers for atmospheric XCO 2 and other trace gases measurements [71]. We compared GM-XCO 2 retrievals to these ground-based measurements to identify potential measurement biases. XCO 2 retrievals from TCCON has a high accuracy with approximately 0.25%, which can be used for validation of satellites observed XCO 2 [72]. In this study, we selected atmospheric CO 2 measurement from TCCON sites with data between 2003.01 to 2016.03 for GM-XCO 2 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 2 comparison. As a result, there were 3603 pairs of GM-XCO 2 and TCCON XCO 2 , from 15 TCCON sites, used for comparison (GM-XCO 2 minus TCCON XCO 2 ). R 2 , mean bias, mean absolute bias, and standard deviation between them were 0.96, 0.20, 0.90, and 1.16, respectively.

Climate Factors
Land surface temperature (LST) and precipitation were used for evaluating the effects from climate change in XCO 2 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 2 .

Trends in Regional XCO2 Summer Drawdown
Here we identify the trends in GM-XCO2 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.

Satellite Observed XCO2 Trends after Removing Global Mean Growth
As a first attempt, we removed the NOAA global mean XCO2 growth rate from the GM-XCO2 in order to calculate summer drawdown trends using the mean value of detrended GM-XCO2 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 −1 for decreasing grid cells). Blue grid cells with black outlines, indicating significant increasing trends, 0.05 ± 0.05 ppm yr −1 , were observed in some urban regions like Beijing and Shanghai.
Regions around cites like Beijing and Shanghai showed increasing XCO2 at more extreme rates than the global mean. Blue grid cells indicate CO2 emission increases (0.045 ± 0.055 kg C m −2 yr −2 ) near urban centers of Eastern Asia in Figure 5a. This is because fossil emissions have increased faster in

Trends in Regional XCO 2 Summer Drawdown
Here we identify the trends in GM-XCO 2 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.

Satellite Observed XCO 2 Trends after Removing Global Mean Growth
As a first attempt, we removed the NOAA global mean XCO 2 growth rate from the GM-XCO 2 in order to calculate summer drawdown trends using the mean value of detrended GM-XCO 2 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 −1 for decreasing grid cells). Blue grid cells with black outlines, indicating significant increasing trends, 0.05 ± 0.05 ppm yr −1 , were observed in some urban regions like Beijing and Shanghai. 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 CO2 emission domes in order to isolate trends in local biospheric influence on XCO2.
(a) XCO2 trends in GP corrected with NOAA (b) Fossil fuel CO2 emission trend from ODIAC

Satellite-Observed XCO2 Trends after Removing Influence of the Urban CO2 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-XCO2 observations: GM-XCO2 fitting (Figure 6a), CT-XCO2 fitting (Figure 6b), and GEOS-Chem simulated XCO2 fitting at two different spatial scales: a coarse 2 × 2.5 degree lat/lon simulation and a highresolution nested simulation for Eastern Asia (Figure 6c-d). All methods showed GM-XCO2 residuals decreasing in many EC grid cells from 2003 to 2015 indicating greater CO2 uptake during the GP. For the entire EC domain, detrending the XCO2 time series of each grid by the long-term trend of each grid cell (GM-XCO2 detrending using the full record including growing season and dormant season) resulted in a small trend in the corrected GM-XCO2 during the GP, −0.058 ± 0.029 ppm yr −1 . 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 XCO2 from remote and local non-biospheric CO2 fluxes resulted in larger GP CO2 uptake trends: −0.084 ± 0.090, −0.081 ± 0.088, and −0.070 ± 0.093 ppm yr −1 (mean ± 1 standard deviation in space), for CT-XCO2 and GEOS-Chem coarse and fine resolution, respectively. Of all the grid cells, larger GP CO2 uptake trends were observed in 80%, 79%, and 77% of the grid cells depending on the correction method applied, CT-XCO2 fitting, and GEOS-Chem fittings at coarse and fine resolution, respectively. Grid cells with statistically significant (p < 0.05) increasing GP CO2 uptake trends occupied 40%, 41%, and 39% of EC, respectively. The trends for just the grid cells that showed increased CO2 uptake were −0.116 ± 0.067, −0.113 ± 0.069, and −0.108 ± 0.068 ppm yr −1 , 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 CO2 uptake in EC due to biological activity. We consider the corrected GP trend using the GM-XCO2 annual increase of −0.058 ppm yr −1 to be the lower estimate of biospheric Regions around cites like Beijing and Shanghai showed increasing XCO 2 at more extreme rates than the global mean. Blue grid cells indicate CO 2 emission increases (0.045 ± 0.055 kg C m −2 yr −2 ) 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 2 emission domes in order to isolate trends in local biospheric influence on XCO 2 .

Satellite-Observed XCO 2 Trends after Removing Influence of the Urban CO 2 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 2 observations: GM-XCO 2 fitting (Figure 6a), CT-XCO 2 fitting (Figure 6b), and GEOS-Chem simulated XCO 2 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 2 residuals decreasing in many EC grid cells from 2003 to 2015 indicating greater CO 2 uptake during the GP. For the entire EC domain, detrending the XCO 2 time series of each grid by the long-term trend of each grid cell (GM-XCO 2 detrending using the full record including growing season and dormant season) resulted in a small trend in the corrected GM-XCO 2 during the GP, −0.058 ± 0.029 ppm yr −1 . 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 2 from remote and local non-biospheric CO 2 fluxes resulted in larger GP CO 2 uptake trends: −0.084 ± 0.090, −0.081 ± 0.088, and −0.070 ± 0.093 ppm yr −1 (mean ± 1 standard deviation in space), for CT-XCO 2 and GEOS-Chem coarse and fine resolution, respectively. Of all the grid cells, larger GP CO 2 uptake trends were observed in 80%, 79%, and 77% of the grid cells depending on the correction method applied, CT-XCO 2 fitting, and GEOS-Chem fittings at coarse and fine resolution, respectively. Grid cells with statistically significant (p < 0.05) increasing GP CO 2 uptake trends occupied 40%, 41%, and 39% of EC, respectively. The trends for just the grid cells that showed increased CO 2 uptake were −0.116 ± 0.067, −0.113 ± 0.069, and −0.108 ± 0.068 ppm yr −1 , respectively (mean ± 1 standard deviation in space).

Trends in Flux Estimates
The corrected GM-XCO2 trends of −0.070 to −0.084 ppm yr −1 in EC suggest an increased biospheric uptake. In order to estimate the magnitude of the surface CO2 fluxes that would produce that XCO2 anomaly, we used the CT and MACC models to relate CO2 fluxes to atmospheric CO2 concentration anomalies. We explored correlation between model simulated XCO2 changes daily or over 16 days and net CO2 flux changes in the global or EC region shown in The similar pattern of trends of these last three correction methods showed that there was robust evidence of increased GP CO 2 uptake in EC due to biological activity. We consider the corrected GP trend using the GM-XCO 2 annual increase of −0.058 ppm yr −1 to be the lower estimate of biospheric influence on EC CO 2 uptake because the local biospheric influence is in the GM-XCO 2 signal itself and a portion of it is removed by the fitting. CT-XCO 2 should depict the general pattern of atmospheric CO 2 concentration over city areas, as the posterior CO 2 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 2 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.

Trends in Flux Estimates
The corrected GM-XCO 2 trends of −0.070 to −0.084 ppm yr −1 in EC suggest an increased biospheric uptake. In order to estimate the magnitude of the surface CO 2 fluxes that would produce that XCO 2 anomaly, we used the CT and MACC models to relate CO 2 fluxes to atmospheric CO 2 concentration anomalies. We explored correlation between model simulated XCO 2 changes daily or over 16 days and net CO 2 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 2 = 0.84, p < 0.01) and MACC (Figure 7d, R 2 = 0.49, p < 0.01) with slopes approaching 0.26 ppm/(g C/m 2 /day) from CT and 0.21 ppm/(g C/m 2 /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 2 is conserved during atmospheric transport within the global domain, whereas at the regional domain, transport can be a source or sink of CO 2 . 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 2 and vegetation parameters) became stronger both for CT (Figure 7c from CT and 0.21 ppm/(g C/m 2 /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, CO2 is conserved during atmospheric transport within the global domain, whereas at the regional domain, transport can be a source or sink of CO2. We estimated the CO2 flux change from the corrected GM-XCO2 trends and the model correlation between XCO2 16-day changes and CO2 flux change shown in Table 2. The total XCO2 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-XCO2 fitting, background GEOS-XCO2 fitting, and nested GEOS-XCO2 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 2 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 2 using the MACC transform ratio. That is basically consistent with GPP increase of 9.24 g C m −2 We estimated the CO 2 flux change from the corrected GM-XCO 2 trends and the model correlation between XCO 2 16-day changes and CO 2 flux change shown in Table 2. The total XCO 2 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 2 fitting, background GEOS-XCO 2 fitting, and nested GEOS-XCO 2 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 2 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 2 using the MACC transform ratio. That is basically consistent with GPP increase of 9.24 g C m −2 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 12 m 2 shown in Figure S8. Table 2. GM-XCO 2 decreasing trends and estimated CO 2 uptake changes from 2003 to 2015 during the vegetation growing period over Eastern China (EC). There are results from three different XCO 2 background corrections and two flux-XCO 2 transform ratios. XCO 2 background corrections are shown in Section 3.2. Transform ratios are shown in Figure 7. Uncertainties in XCO 2 trends and changes refer to 1 standard deviation of spatial grid-scale variability in EC. 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 2 flux during the GP during 2003 to 2016 over EC was −99.8 ± 30.6 g C m −2 yr −1 (uncertainty is 1 standard deviation of inter-annual variability). Annual fossil fuel CO 2 emission for the same period was 328.48 ± 91.38 g C m −2 yr −1 (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.

Possible Contribution from Climate Variability, Independent of Land Management
To check the relationship between XCO 2 and climate change, we present the Pearson correlation coefficients between corrected XCO 2 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 2 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 2 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 2 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 2 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].

Discussion
Vegetation greening trends ( Figure 4) and XCO2 decreasing trends ( Figure 6) during the GP show similar patterns of enhanced greening and CO2 uptake in Eastern China (EC). This indicates that enhancements in vegetation activity have had a measurable influence on regional atmospheric XCO2 concentrations due to enhanced CO2 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 XCO2 trend, and consider the uncertainty in our analysis.

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].

Uncertainty of GM-XCO2 Trends
The GM-XCO2 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 XCO2 in each grid cell over EC ( Figure S11

Discussion
Vegetation greening trends ( Figure 4) and XCO 2 decreasing trends ( Figure 6) during the GP show similar patterns of enhanced greening and CO 2 uptake in Eastern China (EC). This indicates that enhancements in vegetation activity have had a measurable influence on regional atmospheric XCO 2 concentrations due to enhanced CO 2 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 2 trend, and consider the uncertainty in our analysis.

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].

Uncertainty of GM-XCO 2 Trends
The GM-XCO 2 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 2 in each grid cell over EC ( Figure S11 Here, we show the temporal change of integrated XCO 2 (the merging of data from the three sensors) and GM-XCO 2 (gap-filled) using available data in EC as shown in Figure 9. GM-XCO 2 agreed well with integrated XCO 2 with a mean difference of 0.25 ± 1.33 ppm. Some differences exist during summer for distinct XCO 2 spatial variations shown with mean ± 1 STD GM-XCO 2 (red dashed lines). We found no significant trend in the GP residuals between integrated XCO 2 and the gap-filled GM-XCO 2 using the Mann-Kendal trend test. Therefore, we conclude that the XCO 2 decreasing trends are not caused by the gap-filling method. sensors) and GM-XCO2 (gap-filled) using available data in EC as shown in Figure 9. GM-XCO2 agreed well with integrated XCO2 with a mean difference of 0.25 ± 1.33 ppm. Some differences exist during summer for distinct XCO2 spatial variations shown with mean ± 1 STD GM-XCO2 (red dashed lines). We found no significant trend in the GP residuals between integrated XCO2 and the gap-filled GM-XCO2 using the Mann-Kendal trend test. Therefore, we conclude that the XCO2 decreasing trends are not caused by the gap-filling method.
In addition, we checked the GM-XCO2 product against TCCON measurements as shown in Figure S5. The mean absolute difference (bias) between GM-XCO2 and TCCON observations with 3608 co-located data pairs was 0.90 ± 1.16 ppm. Individual grid-cell GM-XCO2 bias followed a Gaussian distribution, therefore the XCO2 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 XCO2 total change (0.840-1.008 ppm from 2003 to 2016). XCO2 differences between GM-XCO2 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-XCO2 trends over time.   Red line is the regional mean (solid) GM-XCO 2 ± 1 STD (dashed). Dark green points are the difference between mean GM-XCO 2 and median integrated XCO 2 (+355 ppm).
In addition, we checked the GM-XCO 2 product against TCCON measurements as shown in Figure S5. The mean absolute difference (bias) between GM-XCO 2 and TCCON observations with 3608 co-located data pairs was 0.90 ± 1.16 ppm. Individual grid-cell GM-XCO 2 bias followed a Gaussian distribution, therefore the XCO 2 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 2 total change (0.840-1.008 ppm from 2003 to 2016). XCO 2 differences between GM-XCO 2 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 2 trends over time. Here, we show the temporal change of integrated XCO2 (the merging of data from the three sensors) and GM-XCO2 (gap-filled) using available data in EC as shown in Figure 9. GM-XCO2 agreed well with integrated XCO2 with a mean difference of 0.25 ± 1.33 ppm. Some differences exist during summer for distinct XCO2 spatial variations shown with mean ± 1 STD GM-XCO2 (red dashed lines). We found no significant trend in the GP residuals between integrated XCO2 and the gap-filled GM-XCO2 using the Mann-Kendal trend test. Therefore, we conclude that the XCO2 decreasing trends are not caused by the gap-filling method.
In addition, we checked the GM-XCO2 product against TCCON measurements as shown in Figure S5. The mean absolute difference (bias) between GM-XCO2 and TCCON observations with 3608 co-located data pairs was 0.90 ± 1.16 ppm. Individual grid-cell GM-XCO2 bias followed a Gaussian distribution, therefore the XCO2 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 XCO2 total change (0.840-1.008 ppm from 2003 to 2016). XCO2 differences between GM-XCO2 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-XCO2 trends over time.

Uncertainties in Constraining the Urban Dome and Remote Influence on XCO 2 Trends
The XCO 2 trends observed over EC result from a combination of remote and local influences. Accurately quantifying the influence of regional biospheric CO 2 fluxes in EC depends on how well we remove the influence from non-local biospheric CO 2 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 2 increase to better study the local biospheric component of XCO 2 changes. The first method using the NOAA-reported global mean CO 2 is not able to remove the influence from the local urban CO 2 dome, which increases over time and therefore underestimated the trend in XCO 2 due to local biospheric influence. In addition, the time series fitting of GM-XCO 2 removed the local annual biospheric contribution to the XCO 2 trend and underestimated the signal of interest, although there was still a component of GS trends in XCO 2 that served as a lower bound of local biospheric influence. The trends calculated from detrending using simulated XCO 2 using CT and GEOS-Chem transport (with and without nested grid cells) of CT surface CO 2 fluxes were likely the best representation of local biospheric influence on XCO 2 because they are designed to remove trends caused by the remote biospheric CO 2 fluxes over time, but retain trends in local biospheric CO 2 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 2 dome was not critical to defining GP trends. Therefore, XCO 2 decreasing trends of −0.084 ± 0.090, −0.081 ± 0.088, and −0.070 ± 0.093 ppm yr −1 were our best estimates of XCO 2 change caused by EC region local biospheric CO 2 fluxes with uncertainty listed is the grid-scale spatial variability in trends ( Table 2). These trends result in a total XCO 2 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 2 fluxes could result in artificial trends in local biospheric CO 2 fluxes in the study region, however it is very unlikely that such an artificial signal from remote CO 2 flux uncertainty would correspond to the exact region that shows vegetation greening. Second, an overestimate in local fossil fuel CO 2 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 2 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 2 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 2 emissions uncertainty. However, the magnitude of the inferred sink is sensitive to local CO 2 emissions uncertainty. This demonstrates the importance of independent constraints of fossil and biospheric CO 2 fluxes in order to understand CO 2 flux trends near urban centers.

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 −2 over EC. To understand the significance of these changes, we compare with mean posterior flux estimates from CT ( Table 3). 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 2 emissions. By comparison, the increased carbon uptake during the GP was slightly underestimated in the interannually varying CT posterior CO 2 fluxes, which is about 6.4% to 10.3% of regional CT posterior annual biopsheric flux.
In this estimation, we calculated the biospheric CO 2 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 2 trends in the dormant season, after CT-XCO 2 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. Table 3. Comparison of the estimated carbon uptake increase (6.36 to 10.30 g C m −2 ) due to greening with various posterior CO 2 fluxes inferred from CT 2017 for the region during the growing period (GP) and the whole year. Posterior CO 2 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 2 sink enhancement as percentages of the various GP CT CO 2 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 ) 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 2 fluxes to atmospheric CO 2 observations. Our approach may likely underestimate the surface CO 2 sink strength as some of the negative CO 2 anomaly is transported out of the study domain. However, it is unlikely that an anomaly caused by a CO 2 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.

Conclusions
Significant greening trends in Eastern China were observed by various parameters including NDVI trends of 0.70 ± 0.61% yr −1 (MODIS) and 0.71 ± 0.75% yr −1 (AVHRR) and an EVI trend of 0.86 ± 0.60% yr −1 (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 2 uptake in this region.
Satellite-observed XCO 2 trends showed a corresponding increase in the carbon sink in the region of greening in EC. Significant XCO 2 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 2 trends across EC varied from −0.070 to −0.084 ppm yr −1 . Estimated local CO 2 uptake increases were 6.36 to 10.30 g C m −2 over EC during the growing period, which is consistent with a MODIS observed GPP increase of 9.24 g C m −2 . 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 2 flux change.