VPRM-CHINA : Using the Vegetation , Photosynthesis , and Respiration Model to partition contributions to CO 2 measurements in Northern China during the 2005-2009 growing seasons

Accurately quantifying the spatiotemporal distribution of the biological component of CO2 surface-atmosphere exchange is necessary to improve top-down constraints on China's anthropogenic CO2 emissions. We provide hourly fluxes 15 of CO2 as Net Ecosystem Exchange (NEE; μmolCO2ms) on a 0.25°x0.25° grid by adapting the Vegetation, Photosynthesis, and Respiration Model (VPRM) to the eastern half of China for the time period from 2005-2009; the minimal empirical parameterization of the VPRM-CHINA makes it well-suited for inverse modeling approaches. This study diverges from previous VPRM applications in that it is applied at large scale to China’s ecosystems for the first time, incorporating a novel processing framework not previously applied to existing VPRM versions. In addition, the VPRM20 CHINA model prescribes methods for addressing dual-cropping regions that have two separate growing season modes applied to the same model gridcell. We evaluate the VPRM-CHINA performance during the growing season and compare to other biospheric models. We calibrate the VPRM-CHINA with ChinaFlux and FluxNet data and scale up regionally using Weather Research and Forecasting (WRF) Model v3.6.1 meteorology and MODIS surface reflectances. When combined with an anthropogenic emissions model in a Lagrangian particle transport framework, we compare the ability of VPRM25 CHINA relative to an ensemble mean of global hourly flux models (NASA CMS) to reproduce observations made at a site in Northern China. The measurements are heavily influenced by the Northern China administrative region. Modeled hourly timeseries using vegetation fluxes prescribed by VPRM-CHINA exhibit low bias relative to measurements during the MaySeptember growing season. Compared to NASA CMS subset over the study region, VPRM-CHINA agrees significantly better with measurements. NASA CMS consistently underestimates regional uptake in the growing season. We find that 30 during the peak growing season, when the heavily cropped North China Plain significantly influences measurements, VPRM-CHINA models an CO2 uptake signal comparable in magnitude to the modeled anthropogenic signal. In addition to demonstrating efficacy as a low-bias prior for top-down CO2 inventory optimization studies using ground-based 5 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504 Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c © Author(s) 2017. CC BY 4.0 License.


Introduction
In 2006, China surpassed the USA as the world's leading anthropogenic carbon dioxide (CO 2 ) emitter.China's contribution to world CO 2 emissions has been growing steadily and now constitutes approximately 26% of the world total, compared to the USA's 17%, accounting for 60% of the overall growth in global CO 2 emissions over the past 15 years (EIA, 2017).China and the USA made an historic joint announcement on national carbon commitments in November 2014, an unprecedented 10 form of political coordination by the two countries to advance United Nations climate negotiations.In addition, China pledged at the 2015 UN climate summit in Paris to peak carbon emissions by 2030; in March 2016, China released its 13 th Five-Year Plan to strengthen its strategies to achieve its emission targets (Tollefson, 2016).An accurate assessment of CO 2 fluxes within China is not only a critical advance in testing the bottom-up emission inventories that provide the baselines for setting such policy commitments and measuring progress, it also broadens understanding of the country's contributions to 15 climate change beyond the sources conventionally targeted for control.Eventually such observation-based assessments might be formally integrated into regulatory processes to strengthen baselines, to widen the scope of control, and to assess policy progress and compliance.
Good prior estimates of the spatiotemporal structure of CO 2 surface exchanges are needed to reduce the uncertainty in top-20 down optimizations where atmospheric observations are used as a constraint to improve bottom-up flux inventories.As a first step towards evaluating China's anthropogenic emission inventories on an intra-annual basis, it is necessary to also model vegetation contributions to the CO 2 signal during the growing season.Previous studies (Wang et al., 2010) relied on CO 2 to CO ratios to estimate annual anthropogenic CO 2 emission enhancements from winter-time observational data alone.
The large diurnal CO 2 uptake and emission vegetation signal in the growing season complicated modeling the anthropogenic 25 CO 2 signal during these times of the year.Seasonal variations in anthropogenic emissions patterns from both shifts in atmospheric transport and emission sources themselves are therefore insufficiently captured when dormant season observations alone are used to estimate annual emissions.
This study adapts the Vegetation, Photosynthesis, and Respiration Model (VPRM; Mahadevan et al., 2008) to model CO 2 30 vegetation fluxes on an hourly basis on a 0.25°x0.25°grid from 2005-2009.VPRM-CHINA is empirically driven with very low dimensional parameterizations, making it particularly suitable as a biogenic CO 2 flux prior model for top-down inverse Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.
3 model analysis of China's emissions.We demonstrate its validity as a prior for use in atmospheric inversions that constrain anthropogenic CO 2 emissions on inter-and intra-annual scales by comparison to hourly observations at a site in Miyun, China, 100km northeast of the Beijing urban center.

Methods 5
As described in detail by Mahadevan et al. (2008), modeled CO 2 vegetation fluxes (Net Ecosystem Exchange; NEE) from VPRM-CHINA are calibrated against observed CO 2 fluxes for each dominant vegetation class represented in the study domain.CO 2 exchange is dependent on ecosystem temperature and sunlight, driven here by high-resolution meteorology from the Weather Research and Forecasting (WRFv3.6.1)model (http://wrf-model.org).The photosynthetic capacity of ecosystems is controlled by vegetation greenness and water availability, and these factors are obtained from MODIS remote-10 sensing (https://lpdaac.usgs.gov)land cover and surface reflectance data sets.We define uptake and release of carbon relative to the atmosphere such that the photosynthesis (Gross Primary Productivity; GPP) term is negative (representing CO 2 uptake from the atmosphere) and the ecosystem respiration (R eco ) term is positive (representing CO 2 release to the atmosphere).Modeled CO 2 in the growing season is ultimately evaluated against hourly averaged observations collected from 2005 to 2009. 15 This study diverges from previous VPRM applications in three main ways: (1) the VPRM-CHINA is applied at large scale to China's ecosystems for the first time; (2) the scaling involves a novel processing framework not previously applied to existing versions of the VPRM; and (3) the VPRM-CHINA model prescribes methods for addressing winter wheat and corn dual-cropping regions that have two growing season modes for the same pixels.
We begin this section with an overview of the observational record used in this study (Sect.2.1).Sect.2.2 presents details of VPRM-CHINA, including processing of model drivers and model calibration.Sect.2.3 introduces the anthropogenic emissions inventory used in this study.Sect.2.4 concludes the methods section, summarizing the derivation of the modeled 25 CO 2 time-series.

Model Overview
We follow the general model framework established by Mahadevan et al. (2008)

×𝐸𝑉𝐼×𝑃𝐴𝑅 + 𝛼×𝑇 + 𝛽
(1) The parameters l, a, b and PAR 0 are empirically adjusted based on calibrations against observed NEE from eddy flux data for each MODIS vegetation class based on the International Geosphere-Biosphere Programme (IGBP; MCD12Q1) in the 15 domain.PAR 0 represents the half-saturation value of photosynthetically active radiation.In addition, we set a minimum temperature threshold for each vegetation class, T=T low (1≤T low ≤5°C) prescribing a baseline of soil respiration at very low air temperatures (Hilton et al., 2013;Mahadevan et al. 2008).T low is derived from fits to site-level data; see Sect.2.2.4 for details.

20
The temperature sensitivity is defined as below, where T min , T max , and T opt represent minimum, maximum, and optimal temperatures for photosynthesis respectively and are set at literature values for each vegetation class.Temperature T is the hourly averaged 10-min surface temperature output from the WRFv3.6.1 meteorological model (Sect.2.2.3).With the exception of winter wheat, we use the same T min , T max , and T opt for each ecosystem type in our domain as in the Mahadevan et al. (2008) North America study.The similarity of latitudes for each ecosystem type between the Mahadevan et al. (2008) 25 study and this study makes this an appropriate approximation.The only ecosystem category in our domain not represented by Mahadevan et al. (2008) is winter wheat; as such we set our winter wheat T min = 0C and T opt = 20C , the lower values relative to other crop types reflecting the lower temperatures of the winter wheat growing season (Acevedo et al. 2002).

15
P scale is set to 0 for water, snow and ice, and unclassified pixels at all times.For evergreen classes at all times and other vegetation classes at maximum greenness we set P scale to 1.We represent phenology as a fraction of LSWI for non-evergreen vegetation classes from (1) onset greenness increase to greenness maximum, and (2) onset greenness decrease to greenness minimum: Ecosystem timing dates for selection of the appropriate P scale parameterization for each pixel were obtained from the MOD12Q2 phenology product, detailed in Sect.2.2.2.

Satellite Data Processing
We use tiles from three MODIS products on a 500m sinusoidal grid to model GPP: 8-day average MOD09A1 surface 25 reflectance bands 2, 6, 1 and 3 representing the near IR, short wave IR, red, and blue regions respectively; annual MCD12Q1 land use categories based on IGBP land classifications; and annual MOD12Q2 ecosystem timing dates.We do not include MODIS surface reflectance data from the Aqua satellite due to failure of a majority of band 6 detectors after launch, which affected availability of high-quality data during the time period investigated in this study.All datasets were downloaded using the Reverb tool in NASA's Earth Observing System Data and Information System (http://reverb.echo.nasa.gov/).We 30 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.
then used MRT to (1) mosaic tiles associated with our spatial domain and (2) reproject to a WGS84 datum Geographic Coordinate system on a 500m grid.
The dominant IGBP ecosystem types represented in the domain are effectively constant over the 5-year study period.
Together, water and the major photosynthesizing land classes constitute 98% of the domain (Figure 1).We follow the 5   2016).We designate all cropland classes south of 32N as rice; all cropland classes between 32N and 38N as winter wheat/corn dual croplands, and all croplands north of 38N as corn.We further justify the designation of the 32N-38N 5 cropland classes as winter wheat/corn by noting a distinct bimodal pattern of average EVI for that region (Fig. 2), similar to Yan et al. (2009).For the dual croplands only, we manually set for each of the two crop modes a multi-year (2005)(2006)(2007)(2008)(2009) mean of phase timings for each degree latitude from 32N to 38N, obtained from EVI averaged across all cropland pixels at the respective degree latitude.10

WRF Temperature and Radiation Fields
Hourly averaged surface temperature (T2) and downward shortwave radiation (SWDOWN) fields were derived from 10minute WRF output.The WRF model was initialized with NCEP FNL 1 °x1 ° resolution reanalysis data (NCEP, 2000) and run for independent 24h periods with three domains spanning the study temporal domain, excluding a 6h spinup time according to practice established by Nehrkorn et al. (2010) and Hegarty et al. (2013).We employ continuous nudging in the 15 outer domain only and we do not nudge any fields within the planetary boundary layer (PBL).The Yonsei University (YSU) PBL scheme is employed.The WRF model is well known to have excess shortwave radiation compared to observations, primarily caused by misrepresentation of clouds and their radiative effects (e.g., Ruiz-Arias et al., 2016).The RRTMG scheme employed in this study includes a method for random cloud overlap in a gridcell (WRF-ARW, 2014).Additional improvements to the treatment of short-wave radiation have been made in more recent versions of WRF than used in this 20 study (Jimenez et al., 2016).Here we will apply an empirical correction to reduce the bias in total incoming shortwave radiation.
The T2 and SWDOWN fields for the outer domain (27km gridcell resolution) are averaged to hourly intervals and then regridded to the same coordinate system as the processed MODIS products using NCL's ESMF bilinear interpolation tool.25 PAR is very closely correlated with shortwave radiation, where SWDOWN ≈ 0.505PAR (Mahadevan et al., 2008).
We quantify the high bias for WRF SWDOWN, and therefore PAR, for our specific study region by comparison of modeled  1).We use unfilled eddy flux data from the Fluxnet 2015 database (http://fluxnet.fluxdata.org/data/fluxnet2015-dataset),ChinaFlux (www.chinaflux.org), and the Fluxnet LaThuile synthesis 5 dataset (http://fluxnet.fluxdata.org/data/la-thuile-dataset/)to calibrate our modeled NEE output for each of our dominant ecosystem types.We used the average of up to nine model pixels of the same IGBP class: NEE from the 500m gridcell nearest to the eddy flux observation site and the surrounding 8 pixels.All observational data were hourly averaged from the original half-hour resolution data sets, and were filtered for sufficiently high frictional velocity, u*, to ensure well-developed turbulence at canopy level prior to calibrating the VPRM-CHINA at each dominant ecosystem type (Goulden et al. 1996).We fit NEE represented in (1) to observations of NEE at each eddy flux site using a non-linear least squares fit (NLS; Gauss-Newton algorithm) in Rv3.2.0.We fit against the subset of non-missing observations and use site-level measurements of air With the exception of dual-cropped croplands, all VPRM-CHINA NEE parameters are obtained from fitting to observations for the whole year.For dual-cropped croplands, we fit the winter wheat and corn modes separately; for rice croplands, used as a validation site only, we use WRF-derived and corrected PAR.With the exception of grasslands, there was insufficient 5 observational data to provide both calibration data and validation data by either a different year in the same site or a different site of the same ecosystem type.In addition, the grassland site used in this study has a history of disturbance.The special cases of grasslands and croplands are discussed further, below.
Within the cropland IGBP ecosystem class we are restricted by availability of observational eddy flux data.Based on USDA 10 agricultural maps we consider corn, winter wheat, and rice as a first order approximation of major croplands influencing CO 2 exchange in the study domain (USDA, 2016).We use 2005 ecosystem data from ChinaFlux site CN-Yuc to calibrate winter wheat/corn dual croplands.Winter wheat and corn parameters were fit to the observational subsets corresponding to times of the year where these crops are prevalent.The winter wheat seasonal subset was defined using dates earlier than July 1, 2005 and the corn seasonal subset was defined using dates on or after July 1, 2005.For corn, the CN-Yuc dataset was manually 15 corrected for a data entry error that began in July of 2005 and lasted through the end of the year.Erroneous instances were flagged where the measured maximum of diurnal NEE uptake lagged measured and modeled PAR by two hours.The NLS fit of NEE was then performed as previously described on this offset-corrected data.Cropland pixels south of 32N were designated as rice and use 2006 ecosystem data from KR-Hae, a rice paddy site in South Korea, to validate rice cropland parameters.KR-Hae dataset contained significant data gaps (43% of data) and could not reliably be used for calibration with 20 the NLS method.Instead we use grassland parameters from Hilton et al. (2013)  transition site (grassland steppe to heavily grazed agricultural) that could impact its ability to represent undisturbed grasslands in the study domain (Zhang et al. 2007;Zhang et al. 2008).

Anthropogenic CO 2 Emissions Inventory
We use the annual anthropogenic CO 2 emissions inventories produced by the Harvard-China Project (ZHAO; Zhao et  season of year changes in activity intensity.We therefore directly scale annual gridcell emissions (as Gg CO 2 ) to µmolCO 2 m -2 s -1 .

WRF-STILT: Derivation of Modeled Hourly CO 2 Timeseries
Each hourly measurement is modeled as advected background air uninfluenced by the study domain plus the integrated effects of CO 2 sources (enhancements relative to background) and sinks (depletion relative to background) from surface 15 processes in the study domain over a specified time period (here, up to seven days back from time of measurement).
We quantify the influence of surface processes using the Stochastic Time-Inverted Lagrangian Transport Model (STILT; Lin et al., 2003), an adjoint that computes surface "footprints" (ppm µmol -1 m -2 s -1 ) for each measurement hour up to 168 hours back from the hour of measurement.An ensemble of 500 hypothetical particles is sent back from the measurement point, 20 driven by WRF meteorology (Nehrkorn et al., 2010).Each footprint ultimately represents the sensitivity of downwind concentration measurements made at a certain hour to upwind surface fluxes.We merge these hourly footprint maps with the vegetation and anthropogenic flux maps pertaining to the appropriate hour (constant in the case of the anthropogenic fluxes).
We can then separately obtain the total anthropogenic enhancement and the vegetation enhancement (or depletion) to the background signal at each measurement hour.25 Background concentrations of CO 2 are estimated using NOAA CarbonTracker CT2015 (NOAA, 2016) provided on a 3hourly 3ºx2º global grid with 25 vertical levels, using a method similar to Karion et al. (2016).A background CO 2 concentration for each hour is calculated as follows: each STILT particle is assigned a background value at the end of its back trajectory.A nearest neighbor approach selects the appropriate CO 2 background concentration based on the particle's 30 end time, latitude, longitude, and altitude.A concentration is only considered to truly represent "background" if at least one Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.12 of the following criteria are met (i) the end point of a particle's back trajectory is the edge of the study's spatial domain; or (ii) if the particle has not reached the edges of the domain, its altitude must be greater than or equal to 3000 masl.A modeled hourly concentration value is then considered valid and included in the analysis if at least 75% of particles satisfy the background selection criteria; the valid background values are then averaged to provide one concentration for each measurement hour.The selection criteria ensure that surface processes in the study domain did not interfere with 5 "background" air during the time period of relevance to the analysis.

3 Results and Discussion
The final VPRM-CHINA product is hourly estimates of NEE on a 0.25°x0.25°grid.Figure 3 displays seasonal averages of hourly VPRM-CHINA NEE over the entire study time period.
The region of high springtime productivity in the North China Plain represents the manually prescribed winter wheat mode.5 Mixed forests at southern latitudes are given the same VPRM-CHINA ecosystem parameters as their only calibration site is a high-latitude mixed forest (CN-Cha; Table 1).This likely leads to an underestimate of mixed forest ecosystem productivity in the south as evidenced by zones of positive summertime mean NEE in southern mixed forest regions.In Sect.3.4, we compare the performance of VPRM-CHINA and CMS in an analysis of multi-year growing season contributions to CO 2 measured at the Miyun station.We conclude with Sect.3.5 where we compare modeled estimates of annual carbon balance within the regions that are estimated as having the greatest influence on the observations.

Seasonal Scaling Factors for Modeled PAR
The WRF PAR bias exhibits seasonal variation, with the highest bias in winter (Fig. 4).We scale modeled PAR in Eq. ( 1) for 10 each season with the derived seasonal scaling factors shown in Fig. 4b, obtained using ranged major axis fits to measured PAR (Legendre and Legendre, 1998).During the dormant season, zones where the P scale phenology term is set to 0 are unaffected by modeled PAR as the light-dependent portion of the NEE equation drops out.We find scaling PAR to be

VPRM-CHINA Calibration Results
Table 2 displays the results from calibration to eddy flux sites in the domain, and compares to results for similar ecosystem 5 classes in North America (Hilton et al., 2013;Mahadevan et al., 2008).Table 3 summarizes the residual standard error (RSE) and the 1-s values from the fits of each parameter.Figure 5   Units are: l: µmolCO 2 m -2 s -1 (µmol PARm -2 s -1 ) -1 ; a: µmolCO 2 m -2 s -1 o C -1 ; b: µmolCO 2 m -2 s -1 ; PAR o : µmol m -2 s -1 5 We set our coordinate system such that uptake by the biosphere is represented as a negative NEE and release to the atmosphere is represented as a positive number.The VPRM-CHINA bias is typically less at hourly timescales (Fig. 5, Table 3), where processes are largely determined by temperature and light parameterizations.However, the unexplained variance is non-random and aggregates over longer timescales (e.g., months and years) and is evident in larger monthly mean biases as in Fig. 6, and in annual sums (Table 4).On annual scales, mean annual temperature and precipitation have been found to 10 dominate carbon exchange in the Asian region (Chen et al., 2013).However, the relationship between cumulative rainfall and LSWI varied depending on whether rainfall is in a high regime or a low regime (Chandrasekar et al., 2010).Future versions of VPRM-CHINA for China would benefit from using Solar Induced Fluorescence (SIF), a more reliable method for quantifying photosynthetic capacity that replaces the EVI, P scale , and W scale terms from the VPRM-CHINA NEE equation (Luus et al., 2017).The time period of this study does not coincide with SIF data availability.At monthly resolutions 15 modeled uptake underestimates observed uptake during the growing season, particularly in the CN-Yuc cropland site during both the winter wheat and corn modes.For all sites and timescales, except evergreen broadleaf and woody savannas that have little influence on the receptor, respiration is underestimated by the model by an additive offset.

Comparison to NASA CMS
We compare VPRM-CHINA performance at hourly timesteps to the NASA CMS weighted ensemble optimal mean of 15 vegetation models (Fisher et al., 2010).The CMS NEE is reported as fluxes of C, provided at 3-hourly resolution on a global 0.5ºx0.5ºgrid.For direct comparison, we convert to fluxes of CO 2 and regrid the CMS using NCLv6.2.1 and a nearestneighbor approach to the same spatial and temporal resolution and extent as the final VPRM-CHINA (hourly; 0.25ºx0.25º).5

Multi-year Growing Season Analysis
Figure 8 displays modeled (ZHAO_VPRM, ZHAO_CMS) and measured CO 2 and residuals for each year during the growing season (May through September) incorporating the peak winter wheat period.We apply VPRM parameters obtained from the previously described calibration.Plots of daily averaged concentrations (Fig. 8a) suggest underestimated uptake by CMS; further examination of modeled-measured residuals at the hourly scale in (Fig. 8b) indicate a systematic underestimate 5 of NEE that is not present in the hourly CO 2 modeled with VPRM-CHINA NEE.Fig. 8c indicates the distributions of ZHAO_VPRM-CHINA CO 2 and measured CO 2 are similar; this is not the case with ZHAO_CMS during times of year where uptake is dominant.
We further examine the relative importance of the vegetation and anthropogenic influence by separately excluding each of 10 the two main components (vegetation and anthropogenic) from the overall unoptimized modeled hourly CO 2 (Fig. 9).

20
We find that accurate modeling of growing season CO 2 in regions of China that are heavily influenced by anthropogenic activity requires incorporation of anthropogenic emissions.While unoptimized CO 2 concentrations modeled only by VPRM were in good agreement with data from aircraft surveys and tower sites in studies in North America (northern New England and Quebec, Matross et al., 2006), the relative magnitudes of the biological and anthropogenic fluxes in the eastern portion of China are of comparable magnitude.In Fig. 9, we compare the ability of VPRM-CHINA and CMS to reproduce growing 5 season CO 2 measurements using a fixed anthropogenic component prescribed by ZHAO.We find that when both vegetation and anthropogenic components are included the VPRM-CHINA vegetation fluxes result in a CO 2 prior that is less biased than one that uses CMS vegetation fluxes.We report mean bias, calculated as the average difference between modeled and measured CO 2 .Based on results displayed in Fig. 8, we conclude that the apparent lower bias of CMS only in Fig. 9d is an artefact of its lower prescription of biosphere uptake in the growing season.10 Extending the results from Fig. 9c and Fig. 9e, we model comparable magnitudes of anthropogenic emissions and biological uptake particularly during the peak growing season.5

Performance of VPRM-CHINA on Annual Timescales
We stress that the VPRM-CHINA is primarily intended as an hourly prior, capturing CO 2 flux covariance with spatial and temporal patterns in vegetation type and status detected by remote sensing.However, analysis at annual timescales is useful to illustrate regional biases.As such, we present analysis of VPRM-CHINA on annual timescales over the study time period by comparison to CMS and Piao et al. (2009).Piao et al. (2009) divided China into nine main regions roughly corresponding 10 to China's administrative regions and examined the vegetation carbon budget on annual timescales using models and vegetation products spanning the time period from 1980 through 2005.
We break down our study region similarly to Piao et al. (2009) as shown in Fig. 11; our study region includes seven of Piao's nine regions.In addition, we calculate the multi-year (2005)(2006)(2007)(2008)(2009) mean annual STILT footprints in the study domain 15 which provides an estimate of the regions that have the greatest influence on the ultimate signal at the Miyun receptor.These  The VPRM-CHINA is poorly constrained by data in Inner Mongolia where we are restricted by the availability of data and use a single degraded grassland site (CN-Du2) to represent over 60% of its landscape.This is particularly problematic as grasslands are shown to have significant ecosystem variation (Zhang et al., 2014).In general, the lack of spatial and, at 25 second order, temporal heterogeneity in both calibration and validation data for the VPRM-CHINA lead to significant error propagation at annual timescales.Jiang et al. (2016) further illustrate this point by evaluating the effect of assimilating more measurements on carbon sink estimates.Jiang et al. (2016) present top-down estimates of land carbon sinks using Carbon Tracker China (CTC) and a Bayesian Inversion (BI) for all of China, noting that the observation network density is highest in the north and east and therefore biased toward constraining exchange in these regions.Increasing the observational 30 constraints from one station to three stations and aircraft data increase carbon sink estimates for the BI and CTC systems by 76% and 95%, respectively (Jiang et al., 2016) and reduce uncertainty in all cases.
In contrast, the VPRM-CHINA is best constrained by data in North and Northeast China.The land categories in these    2014) we find a strong carbon sink in the northeast (Table 1-4), a feature not found by Piao et al. (2009).

Summary and Conclusions
This study shows how the VPRM can be adapted for use in large-scale CO 2 studies in China by (1) streamlining processing 5 of driver datasets at reasonable computational timescales; (2) successfully addressing the significant extent of dual-cropped regions in the North China Plain; (3) calibrating VPRM-CHINA with China-specific observations at the ecosystem level; and (4) demonstrating the low bias of VPRM-CHINA with respect to hourly growing-season CO 2 measured over multiple years at Miyun.We provide hourly estimates of NEE, GPP, and R for the eastern half of China from 2005-2009 on a 0.25°x0.25°grid.We assess its performance in regions that significantly influence CO 2 measured at the Miyun station, and 10 compare it to modeled NEE from other vegetation CO 2 studies across China (NASA CMS, Piao et al., 2009).We use the ZHAO inventory as a control for anthropogenic CO 2 emissions within a WRFv3.6.1-STILTmodeling framework.
The VPRM-CHINA is calibrated with eddy flux data from each major IGBP ecosystem class in the domain (Fig. 1, Table 1).
We separately prescribe a winter wheat mode and corn mode for dual cropland classes within the North China Plain belt 15 (Fig. 2).We find the PAR estimated from WRFv3.6.1 downward shortwave radiation to be overestimated relative to observations by a factor of 1.5 to 2, depending on season (Fig. 4).We find it necessary to scale PAR in Eq. 1 by the seasonal scaling factors for the regional VPRM-CHINA to realistically represent hourly and annual ecosystem carbon fluxes.Overall, the VPRM-CHINA calibration parameters obtained for this study agree with previous studies for similar ecosystem types at similar latitude (Table 2).We find uptake to be well constrained at hourly scales (Fig. 5, Fig. 8, Fig. 9) but underestimated 20 relative to observations at monthly scales (Fig. 6).
We find greater spatial and temporal heterogeneity in VPRM-CHINA relative to CMS over the study time period (Fig. 7).
Accounting for spatial patterns in carbon exchange is necessary for the Lagrangian transport model to capture effects on observed CO 2 at fine grid scales.Indeed, at the hourly scale, the VPRM-CHINA shows significantly greater ability to 25 capture vegetation processes influencing observations at Miyun relative to CMS (Fig. 8, Fig. 9).Due to Miyun's proximity to the North China Plain cropping region, cropland signatures strongly influence partitioning of contributions to modeled CO 2 enhancements relative to CT2015 background concentrations, as evidenced by the strongest drawdown occurring during the peak period of the corn-growing season (Fig. 10).

30
The VPRM-CHINA is well constrained at all timescales by observational ecosystem data in Northeast, North, and South of continuous hourly CO 2 observations (LI-COR Biosciences Li-7000).The site (Miyun; 40°29'N, 116°46.45'E) is in a rural area in Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.Northern China, 100km northeast of the Beijing urban center.Miyun is located south of the foothills of the Yan mountains, and is influenced by clean continental air from the northwest and polluted urban air from the southwest.The vicinity is primarily grasslands, croplands, and mixed temperate forest.Further descriptions of the site and details of the instrumentation of the CO 2 observations are in Wang et al. (2010).

H
, https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.P scale and W scale are functions of the Land Surface Water Index (LSWI).Both LSWI and the Enhanced Vegetation Index (EVI) are derived from the MODIS surface reflectance data set (MOD09A1) as in (3a) and (3b), where the surface reflectance bands used are the red band (r red , band 1); near infrared band (r nir , band 2); blue band (r blue , band 3); and the shortwave infrared band (r swir , band 6).
Mahadevan et al. (2008) convention where (1) we set the NEE of water, urban/built, and snow/ice to zero; and (2) group together (i) savannas and woody savannas; (ii) grasslands, croplands & natural mosaic, and barren & sparse; (iii) deciduous needle-leaf and deciduous broadleaf.The remaining non-dominant vegetation classes not represented in the above (evergreen needle-leaf, closed & open shrublands, permanent wetlands) collectively constitute <1.5% of the total land area and therefore do not appreciably affect the carbon fluxes in the domain.Furthermore, closed and open shrubland ecosystems were found 10 by Mahadevan et al (2008) to be outside of the scope of the VPRM-CHINA due to the inability to adequately capture the influence of inorganic soil carbon pools on observed carbon dioxide fluxes.Pixels corresponding to these ecosystem types have NEE values set to missing.The reflectance data was quality filtered in Rv3.0.2 to accept only the highest quality data under clear sky conditions.15 Inconsistencies in the internal snow-cover flags made it necessary to manually filter erroneous reflectance values due to snow rather than ecosystem photosynthetic activity.

Figure 1 .IIFigure 2 .
Figure 1.Dominant IGBP categories from MOD12Q1 data over study spatial and temporal domain.Circled points represent approximate location and IGBP class of eddy flux sites used in VPRM-CHINA calibration.
PAR to measured PAR at eddy flux sites in the domain and scale our modeled PAR accordingly, by season, using the aggregated PAR across five eddy flux sites.30 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.2.2.4 VPRM-CHINA Calibration Hilton et al. (2013) highlight the importance of tailoring NEE parameter estimates in VPRM to the specific region being studied.As such, we obtain VPRM-CHINA parameters by calibrating the major ecosystem types to representative ecosystem eddy flux data in the domain (Table temperature and PAR.Prior to fitting, we set air temperatures to the respective site T low for instances where the measured air 15 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.temperature falls below that threshold.Initial inputs of l, a, and PAR 0 to the NLS algorithm are based on results from the respective IGBP ecosystem classes in Mahadevan et al. (2008) and Hilton et al. (2013).b is always initially set to 0.
to represent rice cropland parameters and use KR-Hae NEE observations to validate.In addition, KR-Hae did not include PAR observations; therefore we used WRFderived PAR, scaled by seasonal scaling factors in calculations of modeled NEE.Due to the large distances of rice croplands from the Miyun receptor, errors in rice parameterization as a result of this approach are not expected to appreciably affect the final CO 2 concentration estimates at the receptor.25 Grasslands were calibrated using CN-Du2 site data from 2007 and validated using CN-Du2 site data from 2008.The grassland is located in Inner Mongolia and faces low wintertime air temperatures.We find the eddy flux data was collected via an open path system.Insufficient WPL correction used to correct for air density changes (from heating of sensors for example, during the winter time) resulted in possible CO 2 uptake artifacts during the dormant season.Following the 30 convention of Mahadevan et al. (2008) we do not set a T low for grasslands.In addition, the CN-Du2 site represents a Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.
al. 2012) to represent the anthropogenic contribution to the observed CO 2 during the growing season.The ZHAO annual 5 inventories are the first statistically rigorous bottom-up anthropogenic CO 2 inventories for China and integrate data from field studies specific to China's energy processes, technologies, and activity factors with an increased reliance on provinciallevel data relative to national level data.The ZHAO inventories provide emissions in GgCO 2 on a 0.25°x0.25°grid for 2005 and 2009; for 2006-2008 we spatially allocate the total estimated emissions based on the percentage contribution of each gridcell averaged between 2005 and 2009.In addition, we do not apply any temporal activity factors such as time of day or 10

Figure 4 .
Figure 4. Derivation of PAR seasonal scaling factors.(a) Aggregated mean modeled and measured PAR for each eddy flux calibration site by season; error bars are 1-s.(b) Scaling factors for each season, based on fitting modeled to measured PAR using Ranged Major Axis regression.
Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.15 essential to capturing hourly processes influencing measured CO 2 during the growing season (Sect.3.4).Failing to scale WRF-derived PAR with observations would result in a bias factor of 1.5 to 2, leading to an unrealistic overestimate of growing season hourly uptake and annual uptake of CO 2 in the domain.

Figure 5 .
Figure 5. Peak growing season diurnal mean measured NEE (black) along with modeled NEE (red), modeled GPP (purple), and modeled R eco (blue) vs. Local hour of day (UTC+8h).Grey shaded region represents region of positive fluxes (release to atmosphere). .

Figure 6 .
Figure 6.Monthly means of predicted (Modeled) NEE vs Measured NEE at calibration sites, colored by season.Solid line is Standard Major Axis (SMA) regression line; dashed line is 1:1.

Figure 7 Figure 7 .
Figure 7 compares mean and 1-s of annual uptake as kg C m -2 y -1 averaged over the 2005 to 2009 study time period.The CMS product exhibits minimal spatial heterogeneity relative to the VPRM-CHINA product.
Figure 8.Comparison of Modeled and Measured CO 2 during May-September 2005-2009, using two different vegetation models.(a) Daily Averages, for visual clarity; (b) Distribution hourly Modeled-Measured Residuals; (c) Comparison of modeled and measured hourly distributions using a Q-Q plot.

Figure 9 .
Figure 9. Smoothed scatterplots examining relative impact of excluding anthropogenic and vegetation influences in modeling growing season CO 2 .(a) Both hourly anthropogenic and VPRM-CHINA NEE included; (b) both hourly anthropogenic and 3-hourly CMS NEE; (c) hourly VPRM-CHINA NEE only; (d) 3-hourly CMS NEE only; (e) hourly anthropogenic only.All modeled inventories are unoptimized.The blue dashed line represents the 1:1 line; the solid black line represents the SMA fit line described by the equations displayed.Higher point density is represented by darker colors.

Figure 11 .
Figure 11.STILT influences and major administrative regions of China in study domain.Administrative regions are categorized according to convention in Piao et al. (2009).Stippling represents location of 0.25°x0.25°VPRM-CHINA and ZHAO gridcell centers within the regions.Black contour lines display 50 th , 75 th , and 90 th percentiles of mean annual STILT footprints (2005-2009) to highlight regions with the greatest influence on Miyun receptor (solid green circle).Note that two of the regions (Inner Mongolia and Northeast China) are mostly, but not completely, encompassed by our study domain.
regions are appropriately represented by their respective eddy flux calibration sites.For instance, North China has a high 35 prevalence of heavily disturbed grasslands appropriately represented by CN-Du2.In addition, mixed forests in North China are well represented by CN-Cha, which is in close proximity.North and Northeast China fall completely within the 75 th percentile of STILT footprints (Fig.11), thereby increasing confidence in modeled CO 2 during the growing season.Similar toZhang et al. ( China but has sparse representation of its major ecosystem classes in Inner Mongolia, Central, Southeast and Southwest Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-504Manuscript under review for journal Biogeosciences Discussion started: 22 December 2017 c Author(s) 2017.CC BY 4.0 License.This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC.The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE.The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project ofFLUXNET,  20    with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux and AsiaFlux offices.