Analysis of CO2 spatiotemporal variations in China using a weather- biosphere-online-coupled model

Dynamics of atmospheric CO2 has received considerable attention in the literature, yet significant uncertainties 15 remain within the estimates of contribution from terrestrial flux and the influence of atmospheric mixing. In this study we apply the Weather Research and ForecastingWRF-Chem model coupledconfigured with the Vegetation Photosynthesis and Respiration Model (WRF-VPRM) option for biomass fluxes in China to characterize the dynamics of CO2 in the atmosphere. The online coupled WRF-VPRMChem is able to simulate biosphere processes (photosynthetic uptake and ecosystem respiration) and meteorology in one coordinate system. We apply WRF-VPRMChem for a multi-year simulation (2016-2018) 20 with integrated data from a satellite product, flask samplings, and tower measurements to diagnose the spatiotemporal variations of CO2 fluxes and concentrations in China. We find that the spatial distribution of CO2 was dominated by anthropogenic emissions, while its seasonality (with maxima in April 15 ppmv higher than minima in August) was dominated by terrestrial flux and background CO2. Observations and simulations revealed a consistent increasing trend in columnaveraged CO2 (XCO2) of 2.46 ppmv (0.6%/yr) resulting from anthropogenic emission growth and biosphere uptake. WRF25 VPRMChem successfully reproduced ground-based measurements of surface CO2 concentration with mean bias of -0.79 ppmv and satellite derived XCO2 with mean bias of 0.76 ppmv. The model-simulated seasonality was also consistent with observations, with correlation coefficients of 0.90 and 0.89 for ground-based measurements and satellite data, respectively. Tower observations from a background site Lin’an (30.30 ̊N, 119.75 ̊E) revealed a strong correlation (-0.98) between vertical CO2 and temperature gradients, suggesting a significant influence of boundary layer thermal structure on the accumulation and 30 depletion of atmospheric CO2.


Introduction
Climate research requires accurate characterization of atmospheric CO2, which is closely affected by the both atmospheric transport and terrestrial sources and sinks (Bauska et al., 2015;Keenan et al., 2016). Our current knowledge largely comes 35 from interpreting ground-or space-based measurements and model simulations. While observation is limited by spatial and temporal coverages, modelling approaches also suffer from various uncertainties (Shi et al., 2018). Modelling assessment of CO2 is usually conducted through two methods: first, process-or data-driven biosphere models in which terrestrial fluxes are diagnostically calculated with theoretical functions (Tian et al., 2015) or determined through semi-empirical relationships derived from ground measurements and/or satellite products with machine learning techniques (Papale and Valentini, 2003); 40 second, inverse modelling in which prior flux estimates applied in atmospheric transport models are adjusted by observational data and/or satellite products to determine posterior flux (Peylin et al., 2002;Kountouris et al., 2018). Process-driven biosphere models have difficulties capturing spatial and temporal variabilities at fine resolution because parameters calibrated from a limited number of site observations are applied across a variety of land covers (Todd-Brown et al., 2013). Atmospheric inverse modelling is predominantly affected by the presumed prior flux, and different assimilation techniques can give different and 45 even conflicting results (Peylin et al., 2013). These fundamental features highlight the limits of these approaches for accurately modelling carbon dynamics.
Researchers have attempted to reconcile differences between "bottom-up" biosphere models and "top-down" atmospheric inverse models, and recent studies have demonstrated increasing levels of agreement owing to improved understanding of both approaches, such as better parameterization of biosphere processes (Dayalu et al., 2018), more accurately constrained 50 estimates of prior flux Feng et al., 2019), and advanced measurement/satellite instruments that provide high quality data for assimilation (Gaubert et al., 2019); however, critical model disagreements still persist (Kondo et al., 2020).
To bridge the gap between terrestrial flux and atmospheric mixing, a type of weather-biosphere coupled model (Ahmadov et al., 2007;Mahadevan et al., 2008) has been developed to simulate biosphere processes and meteorology conditions in one coordinate system, allowing their interactions to be properly addressed. A few casePrevious modelling studies (Ahmadov et 55 al., 2009;Kretschmer et al., 2012;Park et al., 2018;Beck et al., 2013;Park et al., 2020;Pillai et al., 2012) have demonstrated the potential advantages of coupled weather-biosphere models over pure biosphere/inverse models for short term (a few weeks) simulations, but whether the coupled model is able tocan successfully capture the mesoscale CO2 transport at regional and local scales with significant improvements. But whether it can reproduce the spatial distributions and temporallong-term variations and subsequently estimate carbon fluxes at regional scales with high confidence remains a crucial issue to be 60 addressed.
Understanding the spatiotemporal characteristics of atmospheric CO2 is a key priority in China because of the central role it plays in regulating the climate and environment. In recent years, tremendous efforts have been made in China to control anthropogenic emissions from fossil fuel combustion for both air quality and climate mitigation purposes (Zheng et al., 2018).
While the sources and sinks of air pollutants have been thoroughly examined and well documented (Huang et al., 2020), the 65 dynamics of CO2 at regional to national scales remain poorly understood due to lack of long-term observations and limited modelling studies (Han et al., 2020). Li et al. (2020) applied a weather-biosphere model with tower observations to analyse CO2 fluxes and concentrations over mixed forest and rice paddy in northeast China, but the one-year simulation limits the attempt to investigate interannual CO2 variation which is subject to substantial change (Fu et al., 2019b). Wang et al. (2019) applied satellite products and in-situ observations with inverse modelling to derive posterior carbon fluxes and reported 100% 70 uncertainty for constraining global terrestrial flux. Fu et al. (2020) applied GEOS-Chem simulation with offline Carbon Tracker (Peters et al., 2007) as input to estimate impacts of terrestrial flux and anthropogenic emissions on the annual variation of CO2 concentrations, but regional-scale assessment was limited by coarse grid resolution (2°× 2.5°). Machine-learning technique has also been employed to upscale site observations to regional-scale (Yao et al., 2018;Zhu et al., 2014), but the estimations of carbon budget and dynamics retain large uncertainty due to the diversity of biomass among sites and coarse grid 75 resolution. These pilot studies have shed light on improving the understanding of spatiotemporal characteristics of CO2 in China with modelling or observational methods, but an integrated investigation with both modelling and observations at finescale is urgently needed to expand diagnostic understanding of localized and regional transport, flux, and concentration of CO2 to inform emission management and climate adaption policies (Fu et al., 2019a;Niu et al., 2017;Wang et al., 2019).
In this study we use the Weather Research and ForecastingWRF-Chem model coupledconfigured with the Vegetation 80 Photosynthesis and Respiration Model (WRF-VPRM)Chem) option Mahadevan et al., 2008) to simulate and characterize the spatiotemporal variation of atmospheric CO2 in China from 2016-2018, and also to validate this weatherbiosphere model with recent advanced satellite and tower observations. WRF-VPRMChem has been applied in a few case studies over the United States , Europe , northeast China (Li et al., 2020), and South Korea (Park et al., 2020); this study attempt to apply and evaluate it for a multi-year simulation over China. We first describe 85 the modelling methods and data followed by model validation against observations from multiple datasets, and then present the spatiotemporal variations and estimates of contributions from anthropogenic emissions, terrestrial flux, and background concentrations. Finally, we investigate tower data and reveal the boundary layer thermal structure impacts on atmospheric CO2 accumulation and depletion.

Method 90
We conduct nested WRF-Chem (Version 3.9.1.1)-VPRM) simulations over China (domain shown in Fig.1(a)) and Yangtze River Delta (YRD) region (domain shown in Fig.1(d)) at 20 km and 4 km grid resolution, respectively. Both simulations were configured with 47 vertical layers with model tops at 10hPa. Model configuration in this study followed the work by Hu et al. (2020) and Li et al. (2020). We applied the YSU planetary boundary layer (PBL) scheme (Hong et al., 2006), Morrison microphysics (Morrison et al., 2009), Duhia short-wave radiation (Dudhia, 1989), RRTM long-wave radiation (Mlawer et al., 95 1997), Grell-3 cumulus scheme (Grell and Devenyi, 2002), and Noah land-surface scheme (Chen and Dudhia, 2001), with more details summarized in Table S1. In general, the 4km-grid simulation showed no significant difference as compared to the 20km-grid simulation (demonstrated in Figure S1 and Figure S2), thus the 20km-grid simulation was used to characterize the spatiotemporal distributions of CO2 over China, and the 4km-grid simulation was only used to compare with tower data collected at a background site in YRD. Discussions in the next section will mostly refer to the 20km-grid simulation unless 100 otherwise specified. Initial and lateral boundary conditions for the 20km-grid simulations were derived from the mole fraction product of CarbonTracker (Peters et al., 2007) with 3°×2° resolution. The latest update of column average CO2 (XCO2) concentration assimilation product from CarbonTracker (CT2019) with 1°× 1° resolution (Jacobson, 2020) was also employed to compare with the WRF-VPRMChem simulation. The anthropogenic emission inventory is from the Open-source Data Inventory for Anthropogenic CO2 (ODIAC) with 0.1°× 0.1° resolution (Oda et al., 2018) shown in Fig.1(a); ocean flux 105 is from climatology estimation ). ODIAC has been widely applied in recent modelling studies and demonstrated good agreement with other global inventories (TakahashiHedelius et al., 20092017;Hu et al., 2020). Ocean flux is from climatology estimation (Takahashi et al., 2009); and vegetation fractions and enhanced vegetation index (EVI, shown in Fig.1(b)) are from MODIS (Huete et al., 2002). CO2 from initial and boundary conditions, anthropogenic emission, and terrestrial biogenic flux were tagged as BCG, ANT, and BIO, respectively, to allow the contributions from each process to be identified and quantified 110 through one simulation.
WRF-VPRMChem calculates ecosystem respiration (ER) and gross ecosystem exchange (GEE) with the following functions as: where T is the air temperature at 2m above the surface (T2); , , are vegetation type-dependent parameters; 0 is the vegetation type-dependent half-saturation value of photosynthetically active radiation (PAR); and , , are scaling factors for temperature, water stress, and phenology, respectively. In this study we take the atmosphere as a reference, thus GEE has a negative sign and ER has a positive sign. The current version of WRF-VPRMChem is parameterized ( , , ) for 7 vegetation types ( Fig.1(c)): crops, mixed forest, evergreen forest, deciduous forest, shrub, savanna, and grass. For each 120 modelling grid, ER and GEE are calculated as the weighted averages of each vegetation type based on their fractional abundance. Recent studies Li et al., 2020) have investigated the uncertainty associated with this parameterization through sensitivity simulations and suggested the crops can be further divided into subcategories based on eddy-covariance (EC) tower measurement to improve the model. In this study we used the default parameterization (values presented in Table S2), which has been demonstrated to successfully reproduce the terrestrial flux over northeast China (Li et 125 al., 2020). In contrast, CT2019 applies a process based biosphere model, the Carnegie-Ames Stanford Approach (CASA (Zhou et al., 2020)), driven by year-specific weather and satellite data to simulate terrestrial fluxes (Peters et al., 2007). CASA also estimates photosynthetic uptake based on solar radiation and plant phenology, and estimates respiration as a function of T2.
CASA directly simulates monthly means of Net Primary Production (NPP) and heterotrophic respiration (RH). NPP is the difference between photosynthetic uptake (equivalent to GEE) and autotrophic respiration (RA). The summary of RH and RA 130 is equivalent to ER. Thus, WRF-VPRMChem and CASA are essentially very similar in terms of considering methodology impact; however, it should be noted that to resolve CASA simulated NPP into GEE and RA, CT2019 applies the assumption that GEE is twice that of NPP, which implies that for the same plants the photosynthetic carbon uptake is double the magnitude of autotrophic respiration (but of opposite sign). This assumption is applicable at monthly scale but may have difficulty to reproduce the rapid changes at hourly and daily scales due to impact from weather systems, which will be demonstrated with 135 more details in Section 3.2.
Hourly measurements of CO2 concentrations were collected at the Lin'an Regional Atmospheric Background Station (30.30˚N, 119.75˚E, surroundings shown in Fig.2(a)) with Picarro G1301 and G1302 trace gas analysers mounted on an observation tower at 21 and 55 meters, respectively, above ground level (AGL) and analysed online (data analysis lab shown in Fig.2(b)).
The station is located in the remote area of Hangzhou 138.6 meters above sea level in the middle of a hilly area covered by 140 mixed forest. The observation tower is 60km to the west of downtown center of Hangzhou and 195km to the southwest of Shanghai. Fig.2(c) and (d) presents the wind rose map at Lin'an derived from hourly observations of 10m and 55m wind respectively, which clearly shows the northeast and southwest as prevailing wind directions. The station can properly represent the background atmospheric environment in YRD as demonstrated in previous studies Pu et al., 2020). The tower data provides a representative sampling of CO2 gradients resulting from exchange between atmosphere mixing and 145 terrestrial flux.
Atmospheric samples near the surface were collected at monthly intervals and analysed for CO2 through the National Oceanic and Atmospheric Administration's (NOAA's) Earth System Research Laboratory (ESRL) at four sites (locations shown in December 2016 (Liu, 2018). WRF has been evaluated extensively and consistently performs well for reproducing the meteorology fields and the transport of atmospheric tracers in China (Gao et al., 2015;Tang et al., 2016;Wang et al., 2017;Yang et al., 2019), so this study will only present the simulation performance for CO2 which hasn't been thoroughly discussed in the literature.

Model evaluation
We first evaluate the capability of WRF-VPRMChem to reproduce concentrations of surface CO2 and XCO2, and we find fairly good model performance through the comparison with satellite and ground-based observations. The WRF-VPRMChem simulated surface layer (mid-level height AGL is 12m) CO2 and XCO2 averages between 2016-2018 are demonstrated in Fig.3(a) and (b) respectively. High concentrations were found over industrial areas such as the North China Plain (NCP), Pearl 165 River Delta (PRD), and Yangtze River Delta (YRD), where the surface CO2 and XCO2 were above 440 ppmv and 408 ppmv, respectively; the domain averages were 411 ppmv and 406 ppmv, respectively. While most climate models assume evenly distributed CO2 (Fung et al., 1983;Kiehl and Ramanathan, 1983), our data demonstrates a prominent gradient between industrial and remote areas (e.g., Tibet Plateau, Mongolia), especially for surface CO2, which could be an important source of uncertainty for estimating the long-wave radiation effect (Xie et al., 2018). Spatial patterns of CO2 and XCO2 were in close 170 agreement with ODIAC, indicating the dominant impact of anthropogenic emission in determining the CO2 distribution. WRF-VPRMChem simulated CO2 was generally consistent with CT2019 ( Fig.3(c)), but CT2019 estimated near surface CO2 (midlevel height AGL is 25m) over the coastal industrial areas YRD and PRD because the ocean module used in CT2019 estimated stronger air-sea exchange than the ocean flux by Takahashi et al. (2009) used in WRF-VPRMChem. The two models showed better agreement for XCO2 ( Fig.3 (b) and (e)), but also differed by ~1 ppmv over Taklamakan Desert and along the eastern 175 side of the Tibet Plateau. The OCO-2 swath data were integrated into the corresponding horizontal grids of WRF-VPRMChem and CT2019 respectively, to validate XCO2. Biases of WRF-VPRMChem and CT2019 both fall into the range of ±3 ppmv as shown in Fig.3(c) and (f), respectively, but WRF-VPRMChem apparently provided more details of spatial gradient. WRF-VPRMChem showed well-mixed underestimations and overestimations along neighbouring satellite tracks, while CT2019 tended to overestimate (underestimate) over Tibet Plateau (Taklamakan Desert) where WRF-VPRMChem gave slightly 180 smaller biases. Fig.4(a) and (b) present the raw data pairs between models and OCO-2 with daily interval for WRF-VPRMChem and CT2019, respectively. In general, the WRF-VPRMChem model reproduced OCO-2 well, with mean bias (MB) of 0.76 ppmv, and CT2019 showed MB of 0.54 ppmv, suggesting an overall acceptable performance of the weatherbiosphere model to simulate the spatial distribution pattern of XCO2 in China.
We further analyse WRF-VPRMChem validation against OCO-2 for the seven vegetation types in each season and find no 185 prominent difference (evaluation statistics summarized in Table 1). Regarding vegetation type, the model showed the largest MB over deciduous forest of -1.01 and 1.27 ppmv in summer and winter, respectively, which only covered a very small portion in northeast China. The three most abundant coverage vegetation types in China are grass, crops, and mixed forest. XCO2 simulated by WRF-VPRMChem over grass areas was slightly overestimated by 0.31~0.68 ppmv throughout the year, and the MB over mixed forest was -0.43~0.59 ppmv, indicating a good performance of the model over the vast majority of areas of 190 China. Performance over crops generally showed larger discrepancy than other vegetation types, with MB ranging from 0.66 ppmv in summer to 1.19 ppmv in winter, suggesting the model tends to slightly overestimate column concentration of CO2 over cropland. Li et al. (2020) reported that WRF-VPRMChem underestimated biosphere carbon over rice paddy sites (by ~3%) in northeast China and suggested the parameterization of , , as the most important cause. Cropland differs significantly across China with various types of species such as rice, wheat, and corn, for which literatures reported 195 substantially different rates of ecosystem respiration and photolysis uptake (Gao et al., 2018;Yang et al., 2016;Zhu et al., 2020).
Thus, applying one set of parameters to represent all crops may be responsible for the lingering uncertainty of simulated XCO2.
In terms of seasonal difference, WRF-VRPM showed slightly smaller bias in summer and larger bias in winter, and the correlation coefficients were all ~0.8, consistent with application over the U.S.  which also reported slightly better performance in summer than other seasons. 200 Fig.4 also presents the overall simulation bias against ground-based observations at their raw temporal intervals (monthly for data at ESRL sites, hourly for tower data at Lin'an, and daily for TCCON at Hefei). At the ESRL sites (Fig.4(c)), surface CO2 concentrations were simulated well with minor overestimation by 0.69 ppmv. Evaluation at the Lin'an station was performed with the 4km-grid simulation. The mid-level heights of WRF-VPRM'sChem's first, second, and third layers were 12.3m, 36.9m, and 61.6m, respectively, and simulations were linearly interpolated to 21m and 55m to compare with the tower data. 205 The evaluation at 21m AGL ( Fig.4(d)) shows slight overestimation by 0.02 ppmv, but the evaluation at 55m height ( Fig.4(e)) shows relatively large overestimations by 1.06 ppmv. The discrepancy is likely due to the combined effect of vertical allocation of anthropogenic emission (Brunner et al., 2019) and parameterization of VPRM. Tracer transport models (such as WRF-VPRMChem and CASA) and inverse modelling methods allocate anthropogenic CO2 emission into the near surface layer due to lack of injection height information, which may subsequently lead to systematic overestimation of surface CO2 concentration 210 in industrial areas. Through a regional scale (750×650km) modelling study around the city of Berlin, Brunner et al. (2019) reported that distributing anthropogenic emission into the surface layer overestimated near-surface CO2 concentration by 14% in summer and 43% in winter as compared with considering the vertical profiles of local anthropogenic sources. Lin'an observation tower is located at a densely vegetated area. Validation against OCO-2 suggested that WRF-VPRMChem did not show significantly different performance over different vegetation types as shown in Table 1. As compared to the ESRL 215 background sites which were located in more remote areas with little anthropogenic emission ( Fig.1(a)), Lin'an was more frequently affected by regional anthropogenic emissions which were transported from Shanghai and Hangzhou due to the prevailing northeast wind (Pu et al., 2014), indicating that the emission allocation discrepancy may induce more prominent error at Lin'an. In fact, the 20km-grid WRF-VPRMChem simulation bias at Lin'an were 5.34 and 5.41 ppmv at 21m and at 55m respectively ( Figure S2), significantly larger than the bias at ESRL sites. In addition, both the 20km-grid and 4km-grid 220 simulations showed relatively larger bias at 55m than 21m due to smaller topography roughness and higher wind speed which increases with height according to observations ( Figure S3). CT2019 also substantially overestimated at Lin'an, but the first, second, and third layers' mid-level heights are 25m, 103m, and 247m, respectively, so we do not present the direct comparison with the tower data. Simulated XCO2 from both WRF-VPRMChem and CT2019 were well consistent with the TCCON Hefei site observations as shown in Fig.4(f), with MB by -0.79 ppmv and -0.78 ppmv respectively, and NMB by -0.20% and -0.19% 225 respectively. The 4km-grid simulation showed very similar result to the 20-grid simulation for XCO2 ( Figure S1 and Figure   S2). Recent atmospheric inverse modelling studies (Fu et al., 2019a;Wang et al., 2019;Xie et al., 2018) reported the simulation bias of XCO2 as 0.5-2 ppmv with posterior flux inputs. The WRF-VPRMChem model applied in this study has demonstrated good agreement with the observations though our evaluation.

CO2 seasonal variation and trend in China 230
We next analyse the seasonality of CO2 and XCO2 and find that the terrestrial flux played a more influential role than anthropogenic emission. WRF-VPRMChem successfully reproduced seasonal variations of CO2 at ESRL sites, with a correlation coefficient of 0.90 (Fig.5(a)). The WRF-VPRMChem 4km-grid simulation showed a correlation coefficient of 0.82 with the Lin'an tower observation (averaged for daytime 21m and 55m data). Both the model and measurements showed prominent seasonal cycles for surface CO2 concentrations. The WRF-VPRMChem simulation showed maxima in April (413-235 419 ppmv) and minima in August (399-404 ppmv) as presented in Fig.5(b). The model suggested that the anthropogenic CO2 contribution was 2.6 ppmv in both months, while the biogenic contributions were 3.1 and -1.2 ppmv in April and August, respectively ( Fig.5(d)). Anthropogenic emission (Fig.5(f)) showed a flat curve with relatively higher values in December due to fuel combustion for heating (Zheng et al., 2018). EVI showed maxima in July and August (Fig.5(f)). During summer, photosynthetic uptake almost completely compensated anthropogenic emission, causing the minimum CO2 concentration 240 observed in August, while the higher anthropogenic emission in December and respiration flux in April led to the two corresponding peaks. The anthropogenic XCO2 contributions were 0.5 and 0.6 ppmv in April and August, respectively, and the biogenic contributions were 0.8 and -1.5 ppmv, respectively, suggesting that the seasonality of XCO2 was also primarily dominated by terrestrial flux. Furthermore, the seasonality at high-latitude ESRL sites (UUM and WLG) was stronger than at Lin'an and low-latitude sites (DSI and LLN) because of the larger temperature and photosynthetically active radiation (PAR) 245 gradients. Annual average anthropogenic and biogenic XCO2 contributions were 7.1 and -1.9 ppmv, respectively, indicating that biosphere uptake was an important carbon sink offsetting 27% of anthropogenic emission and slowing the growth of atmospheric CO2.
XCO2 showed similar seasonality, with minima in August and maxima in April and December ( Fig.5(b)). Both WRF-VPRMChem and CT2019 showed good agreement with TCCON Hefei observations with correlations of 0.89 and 0.88, 250 respectively ( Fig.5(e)). However, we note that WRF-VPRMChem simulated drastic changes (e.g., the grey shaded period in Fig.5(e)) that were not reproduced by CT2019. Fig.6 shows the daily concentrations of XCO2 overlaid with horizontal wind speed at 10m AGL from WRF-VPRMChem and CT2019 and highlights large discrepancies over Hefei ( Figure S4 shows the same comparison but using WRF-VPRMChem 4km-grid simulation data). Between April 1 st and 3 rd 2016, an 850 hPa trough associated with a surface cold front moved southeastward from Mongolia to the North China Plain (NCP) (weather maps 255 shown in Fig.6(g)-(i)). At the leading edge of the front, a convergence zone associated with a low pressure center formed, which led to significant cloud formation and subsequently reduced short-wave radiation. As a result, photosynthetic carbon uptake was reduced, leading to enhancement of atmospheric CO2. Meanwhile, the cold front transported anthropogenic CO2 from NCP to YRD, and the convergence zone along YRD ahead of the front facilitated the accumulation of air pollutants and CO2 from anthropogenic emissions. With its coarse spatiotemporal resolution, CT2019 had difficulty reproducing such 260 regional weather systems that can lead to rapid and localized changes in CO2 concentration and terrestrial flux, indicating the importance of fine resolution modelling to better represent the small spatial scale and rapid temporal scale variations of CO2 (Agusti-Panareda et al., 2019).
We also find a notable increasing trend for the 3-year study period. Observed CO2 annual enhancement was 2.2 ppmv/yr (0.56%/yr) at the ESRL sites and 2.3 ppmv/yr (0.54%/yr) at Lin'an. The observed average CO2 concentrations at Lin'an (428 265 ppmv) were substantially higher than those at ESRL sites (407-410 ppmv). The prominent higher levels of CO2 and slightly higher absolute growth rate at Lin'an can be attributed to the influence of the transport regional anthropogenic emission which is growing at rate of 0.82%/yr as suggested by ODIAC. Domain-wide XCO2 was also found to increase by 2.3 ppmv/yr (0.57%/yr) as suggested by OCO-2 and 2.5 ppmv/yr (0.61%/yr) as suggested by the simulation. WRF-VPRMChem reproduced the trends in good agreement with ground and satellite observations. Model simulated budgets suggested that the increasing 270 trends for anthropogenic, biogenic, and background XCO2 were 0.81%/yr, -9.17%/yr, and 0.59%/yr, respectively; the trends for anthropogenic, biogenic, and background CO2 were 4.95%/yr, -0.73%/yr, and 0.59%/yr, respectively. Our findings are consistent with recent measurements and inverse modelling studies but provide process-based estimates for anthropogenic emission and terrestrial flux. Wu et al. (Wu et al., 2012)  although anthropogenic emission increases at a steady rate in East Asia, photosynthetic uptake also serves as an increasing carbon sink due to enhanced EVI (0.29%/yr). However, as the interannual variability (IAV) of terrestrial flux is usually critically large and is affected by both vegetation itself and climate conditions (Fu et al., 2019b;Niu et al., 2017), simulation over longer time periods is necessary in future studies to conclusively comment on the changing trend of biosphere CO2 in 280 China.

Diurnal variation of near-surface CO2 and influence factors
Finally, we examine the diurnal variation of CO2 data at Lin'an station as shown in Fig.7 to reveal the temporal dynamics and atmospheric mixing of CO2 at local scale. While both 21m (Fig.7(a)) and 55m (Fig.7(b)) CO2 show prominent diurnal changes, the variations were larger in summer (JJA) than winter (DJF) and were larger at 21m than 55m, indicating the dominant 285 influence of terrestrial flux over anthropogenic emission in determining the near surface CO2 concentration. Fig.7(c) and (d) present the WRF-VPRMChem simulation bias at 21m and 55m respectively, and Fig.7(e) and (f) present the bias of CT2019 at 21m and 55m respectively. We find that both models prominently overestimated during nighttime, which shall be attributed to the bias in simulating NEE. Li et al. (2020) reported the model overestimated nighttime NEE at a mixed forest site Wuying (47.15˚N, 131.94˚E) by 34% during the growing season (May-Sep.) according to eddy-covariance tower measurement. Fig.7(g)  290 and (h) present the simulated NEE by WRF-VPRMChem and CT2019, respectively, which show close correlations with the CO2 simulation biases. While Lin'an is also covered by mixed forest, our evaluation suggests that WRF-VPRMChem may also overestimate nighttime ecosystem respiration at Lin'an as it has a warmer climate condition than Wuying ( Figure S5), and CT2019 has even greater bias for presenting the diurnal cycles of terrestrial flux.
We also find that planetary boundary layer height (PBLH) significantly affects diurnal accumulation and depletion of 295 atmospheric CO2 as shown in Fig.8(a). During daytime in the growing season, photosynthetic uptake results in lower CO2 concentration; meanwhile, PBLH is also high and allows rapid vertical mixing between near surface and upper air. During nighttime when photosynthesis stops, CO2 from ecosystem respiration starts to accumulate in the shallow stable boundary layer, while the residual layer remains largely decoupled. Thus, atmospheric constituents with surface sources normally exhibit a vertical profile in which concentrations decrease with height in the stable boundary layer Hu et al., 2012). 300 Such boundary layer characteristics are confirmed by CO2 vertical gradients at Lin'an in this study. CO2 at 55m height was consistently lower than the near surface air at 21m during nighttime due to accumulation of respired CO2 in the stable boundary layer. As photosynthetic uptake depleted the near surface CO2 and daytime boundary layer convection developed, the CO2 gradient was gradually weakened from 06:00 to 11:00 LT and remained minimal through the rest of the daytime; at midday when photosynthesis reaches maximum intensity, CO2 at 21m was even lower than at 55m. WRF-VPRMChem roughly 305 reproduced the diurnal profile but noticeably underestimated the intensity of nighttime CO2 difference (ΔCO2) likely due to the bias for simulating night time terrestrial flux as discussed above or underestimation of nighttime boundary layer stability by the PBL scheme (Hu et al., 2012).
The relationship between the near-surface CO2 profile and boundary layer stability is further statistically examined. Fig.8(b) presents the correlation between air temperature gradient (ΔT/ΔH) and CO2 concentration gradient (ΔCO2/ΔH) calculated with 310 diurnal profiles of tower observations averaged for 2016-2018, where ΔT, ΔH, and ΔCO2 is the differences of temperature, height, and CO2 concentration between the two tower levels, respectively. Fig.8(b) clearly demonstrates the influence of boundary layer stability on the CO2 vertical profile, as the correlation between ΔT/ΔH and ΔCO2/ΔH reaches -0.98. On one hand, a more stable PBL with a strongly positive temperature gradient would promote surface CO2 accumulation and lead to a strongly negative CO2 gradient, especially under inversion conditions when upper air has higher temperature (orange area in 315 Fig.8(b)). Conversely, a strongly negative temperature gradient indicates stronger radiation, and subsequently greater photosynthesis and CO2 depletion in the near surface layer, which would result in a positive CO2 gradient (green area in Fig.8(b)) implying a lower CO2 concentration at the surface. While the diurnal variations of ΔCO2 were primarily dictated by local biogenic CO2 fluxes and boundary layer dynamics, the two minor daytime peaks of ΔCO2 at Lin'an, at 10:00 and 18:00 LT ( Fig.8(a)) likely suggest influence of transport of CO2 from urban plumes in the region; for example, from Hangzhou which is 320 60 km away from the tower. Due to rush-hours anthropogenic emissions, CO2 enhancement at Hangzhou relative to a background site exhibited a prominent bimodal curve with two peaks during early morning and early evening (Pu et al., 2018).
Depending on meteorological conditions, particularly wind fields, urban CO2 plumes from cities such as Hangzhou may be transported to the Lin'an site. The influence of boundary layer conditions on CO2 variability has been discussed in several studies through analysis of mountain site ground-based observations (Arrillaga et al., 2019;Esteki et al., 2017;Li et al., 2014), 325 but our study applied tower data as direct evidence to demonstrate the significant impact of PBL thermal structure, which has rarely been documented. More importantly, although WRF-VPRMChem failed to capture the bimodal ΔCO2 peaks at rush hours, because monthly ODIAC data lacked an hourly profile, our analysis reveals the critical importance of careful configuration of the PBL scheme and spatiotemporal distribution of anthropogenic emission for weather-biosphere modelling of atmospheric CO2. 330

Summary and Conclusions
In this study, the spatiotemporal variations of CO2 in China are investigated with measurements from multiple datasets and a weather-biosphere coupled model simulation for 2016-2018. We find consistent higher concentrations over industrial areas with excessive anthropogenic emission and lower concentrations over densely vegetated areas. Observed CO2 concentrations at Lin'an (427 ppmv) are significantly higher than remote ESRL sites (408 ppmv) although they are all established as 335 "background" stations, indicating the dominant influence of anthropogenic emission at regional scales. The Lin'an tower data shows a large negative correlation (-0.98) between vertical CO2 concentration and air temperature gradients, suggesting the significant influence of boundary layer stability on CO2 accumulation and depletion. The online coupled weather-biosphere model WRF-VPRMChem enables process-based estimations of contributions from anthropogenic emission (0.59 ppmv (0.15%)), terrestrial flux (0.16 ppmv (-0.04%)), and background concentration (405.70 ppmv (99.89%)) to average total XCO2. 340 Respective simulation biases of surface CO2 and XCO2 are 0.69 and 0.76 ppmv against ESRL site observations and OCO-2 satellite product with correlations of 0.87 and 0.90, indicating overall good performance of the WRV-VPRMWRF-Chem model. Maximum CO2 concentrations are found in April and minima are found in August for all three years, and the seasonality is reproduced well by the model, which also reveals that terrestrial flux and background concentration dominated the seasonality rather than anthropogenic emission. 345 A steadily increasing trend in XCO2 by 2.46 ppmv (~0.6%/yr) during the study period is demonstrated consistently by both model simulation and satellite product. Budget analysis suggests that anthropogenic emission increased by 0.83%/yr contributing to the 0.81%/yr growth rate of anthropogenic XCO2 enhancement, 27% of which was offset by biosphere uptake.
It is noted that terrestrial flux has significant inter-annual variability, thus a more robust estimation of the terrestrial flux trend should be obtained through a long-term study in the future. The background XCO2, representing contributions from global 350 circulation, increased by 2.37 ppm (0.59%/yr), suggesting that CO2 level in China was growing at the same rate as the rest of the world.
The most significant modelling bias is identified from validation against the Lin'an tower 55m observations, which WRF-VPRMChem 4km-grid simulation overestimated by about 1.06 ppmv with a correlation coefficient of 0.82. The allocation of anthropogenic emission into the surface layer is partially responsible for this modelling bias because Lin'an is closely affected 355 by upwind industrial mega cities in YRD, suggesting the need to include vertical profiles of fossil fuel combustion to properly redistribute the ODIAC for modelling purposes. In addition, diurnal variations of the bias suggest that the modelling discrepancy is also induced by large uncertainty associated with simulating nighttime ecosystem respiration. Representation and parameterization of photosynthetic carbon uptake in VPRM has been continuously improved during the past 10 years since its first release , but ecosystem respiration parameterization is still too simplified to fully represent the 360 autotrophic and heterotrophic respiration of biomass. (Hu et al., 2021). Li et al. (2020) and our study both reveal the urgent need to better calibrate VPRM parameterization over different vegetation types in China, and other methods such as inverse modelling is necessary to further validate the anthropogenic fluxes from ODIAC. Nevertheless, WRF-VPRMChem is demonstrated to be a reliable tool to model the dynamics of CO2 and exchange between the atmosphere and terrestrial flux.
Most importantly, as the online coupled modelling system is able to simulate meteorology and biosphere processes 365 simultaneously, it promotes the opportunity to investigate the interactions between atmospheric mixing and terrestrial flux (Carvalhais et al., 2014;Schimel et al., 2015) while comprehensively considering various factors from both sides that affect CO2 in one coordinate frame, which could be a very helpful tool to support policy makers for balancing short-term carbon cycles at regional scales.

Data availability
The modelling output is accessible by contacting the corresponding author (yjjiang@pku.org.cn, xhu@ou.edu)

Author contributions 375
The concept and ideas to design the integrated simulation and observation analysis are devised by YJ, X-MH, and XD. Simulation is performed by X-MH. OCO-2 satellite product is collected and processed by X-MH. CT2019 assimilation data and ground-based observations are collected by XD. Tower measurement is conducted, processed, and analysed by QM, JP, and YJ. Model evaluation is performed by MY. The manuscript is prepared by XD and X-MH with input and feedback from YJ, MY, QM, JP, and GZ. 380

Competing interests
The authors declare that they have no conflict of interest. , where ̅ is the average of simulations, ̅ is the average of observations. 2 For each season, evaluation statistic with the worst performance (largest absolute value of MB) is highlighted in red, and the one with best performance (smallest absolute value of MB) is highlighted with in bold font. 600