Analysis of CO2 spatiotemporal variations in China using tower data and a weather-biosphere-online-coupled model, WRF-VPRM

Dynamics of CO2 has received considerable attention in the literature, yet significant uncertainties remain within 15 the estimates of contribution from terrestrial flux and the influence of atmospheric mixing. In this study we apply the Weather Research and Forecasting model coupled with Vegetation Photosynthesis and Respiration Model (WRF-VPRM) in China to characterize CO2 dynamics with tower data collected at a background site Lin’an (30.30 ̊N, 119.75 ̊E). The online coupled weather-biosphere WRF-VPRM simulations are able to simulate biosphere processes (photosynthetic uptake and ecosystem respiration) and meteorology in one coordinate system. Simulations are conducted for three years (2016-2018) with fine grid 20 resolution (20 km) to detail the spatiotemporal variations of CO2 fluxes and concentrations. This is the first attempt to apply the weather-biosphere model for a multi-year simulation with integrated data from a satellite product, flask samplings, and tower measurements to diagnose the dynamics of CO2 in China. We find that the spatial distribution of CO2 is determined by anthropogenic emissions, while its seasonality (with maximum concentrations in April 15 ppmv higher than minimums in August) is dominated by terrestrial flux and background CO2. Observations and simulations reveal a consistent increasing 25 trend in column-averaged CO2 (XCO2) of 0.6%/yr resulting from anthropogenic emission growth and biosphere uptake. WRFVPRM successfully reproduces ground-based measurements of surface CO2 concentration with mean bias of -0.79 ppmv (0.20%) and satellite derived XCO2 with mean bias of 0.76 ppmv (0.19%). The model-simulated seasonality is also consistent with observations, with correlation coefficients of 0.90 and 0.89 for ground-based measurements and Orbiting Carbon Observatory-2 (OCO-2) satellite data, respectively. However, evaluation against Lin’an tower data reveals uncertainty within 30 the model for simulating the intensity and diurnal variation of terrestrial flux, which contributes to overestimation by ~5.35 ppmv (1.26%). Lin’an tower observations also reveal a strong correlation (-0.85) between vertical CO2 and temperature https://doi.org/10.5194/acp-2020-1128 Preprint. Discussion started: 2 December 2020 c © Author(s) 2020. CC BY 4.0 License.

Lin'an towerTower observations also revealfrom a background site Lin'an (30.30˚N, 119.75˚E) revealed a strong correlation (-0.8598) between vertical CO2 and temperature gradients, suggesting a significant influence of boundary layer thermal 35 structure on the accumulation and 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 40 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); 45 second, inverse modelling in which prior flux estimates applied in atmospheric transport models are calibrated by observational data and/or satellite products to determine posterior terrestrial flux (Peylin et al., 2002).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 50 modelling is predominantly affected by the presumed prior flux, and different assimilation techniques can give different and 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 55 both approaches, such as better parameterization of biosphere processes (Dayalu et al., 2018), more accurately constrained 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) washas been developed to simulate biosphere processes and meteorology conditions in one 60 coordinate system, allowing their interactions to be properly addressed. A few case studies (Ahmadov et al., 2009;Kretschmer et al., 2012;Park et al., 2018) 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 to reproduce the spatial distributions and temporal variations and subsequently estimate carbon fluxes at regional scales with high confidence remains a crucial issue to be addressed. 65 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 dynamics of CO2 at regional to national scales remain poorly understood due to lack of long-term observations and limited 70 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% uncertainty for constraining global terrestrial flux. Fu et al. (2020) applied GEOS-Chem simulation with offline Carbon 75 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 suffer from coarse grid resolution. These pilot studies have shed light on improving the understanding of spatiotemporal characteristics of 80 CO2 in China with modelling or observational methods, but an integrated investigation with both modelling and observations at fine-scale 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 applyuse the Weather Research and Forecasting model coupled with the Vegetation Photosynthesis and Respiration Model (WRF-VPRM) (Hu et al., 2020;Mahadevan et al., 2008) to simulate and characterize the spatiotemporal 85 variation of atmospheric CO2 in China from 2016-2018, and also to validate this weather-biosphere model with recent advanced satellite and tower observations. WRF-VPRM has been applied in a few case studies over the United States (Hu et al., 2020), Europe (Kretschmer et al., 2012), northeast China (Li et al., 2020), and South Korea (Park et al., 2020); this study is the first attempt to apply and evaluate it for a multi-year simulation at fine scale (20 km) over China. We first describe the modelling methods and data employed followed by model validation against observations from multiple datasets, and then 90 present the spatiotemporal variations and estimates of contributions from anthropogenic emissions, terrestrial flux, and background concentrations. Finally, we probe intoinvestigate tower data and reveal the boundary layer thermal structure impacts on atmospheric CO2 accumulation and depletion.

Method
The WRF-VPRM simulation in this study is configured with 48 vertical layers and 20 km grid resolution. Initial and boundary 95 conditions areWe conduct nested WRF(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 100 (Mlawer et al., 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 105 unless 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 is (Jacobson, 2020) was also employed to compare with the WRF-VPRM 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 is from 110 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 arewere tagged as BCG, ANT, and BIO, respectively, to allow the contributions from each process to be identified and quantified through one simulation.
WRF-VPRM calculates ecosystem respiration (ER) and gross ecosystem exchange (GEE) with the following functions as: 115 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, 120 thus GEE has a negative sign and ER has a positive sign. The current version of WRF-VPRM is parameterized ( , , ) for 7 vegetation types ( Fig.1(c)): crops, mixed forest, evergreen forest, deciduous forest, shrub, savanna, and grass. For each modelling grid, ER and GEE are calculated as the weighted averages of each vegetation type based on their fractional abundance. Recent studies (Hu et al., 2020;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 125 eddy-covariance (EC) tower measurement to improve the model. In this study we applyused the default parameterization, (values presented in Table S2), which has been demonstrated to successfully reproduce the terrestrial flux over northeast China (Li et al., 2020). In contrast, CT2019 applies a pureprocess based biosphere model, the Carnegie-Ames Stanford Approach (CASA ), 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 130 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 is equivalent to ER. Thus, WRF-VPRM 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 135 of autotrophic respiration (but of opposite sign). This assumption is applicable at monthly scale but may contribute to have difficulty reproducingto reproduce the rapid changes at hourly and daily scales due to impact from weather systems, which will be demonstrated with more details in Section 3.2.
MeasurementsHourly measurements of CO2 concentrations arewere collected at the Lin'an Regional Atmospheric Background Station (30.30˚N, 119.75˚E, surroundings shown in Fig.1(d2(a)) with Picarro G1301 and G1302 trace gas analysers mounted 140 on an observation tower at 21 and 55 meters, respectively, above ground level (AGL) and analysed online (data analysis lab shown in Fig.1(e2(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 mixed forest. The hourly Lin'an station tower measurements collected between 2016-2018 provide a representative sampling of theThe 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 145 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 terrestrial flux.  . The TCCON-Hefei site was located in the northwestern rural area of Hefei city and measurements were conducted from September 2015 to December 2016 (Liu, 2018). WRF has been evaluated extensively and consistently performs well for reproducing the meteorology fields and the transport of atmospheric 160 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 only which hasn't been thoroughly discussed in the literature.

Model evaluation
We first evaluate the capability of WRF-VPRM to reproduce concentrations of surface CO2 and XCO2, and we find fairly 165 good model performance through the comparison with satellite and ground-based observations. The WRF-VPRM simulated surface layer (mid-level height AGL is 12m) CO2 and XCO2 averages between 2016-2018 are demonstrated in Fig.23(a) and (b) respectively. High concentrations arewere found over industrial areas such as the North China Plain (NCP), Pearl River Delta (PRD), and Yangtze River Delta (YRD), where the surface CO2 and XCO2 arewere above 440 ppmv and 408 ppmv, respectively; the domain averages arewere 411 ppmv and 406 ppmv, respectively. While most climate models assume evenly 170 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 arewere in close agreement with ODIAC, indicating the dominant impact of anthropogenic emission in determining the CO2 distribution. WRF-VPRM simulated CO2 iswas generally consistent with CT2019 ( Fig.2(d3(c)), but CT2019 estimates lowerestimated near 175 surface CO2 (mid-level height AGL is 25m) over the coastal industrial areas YRD and PRD because the ocean module used in CT2019 estimatesestimated stronger air-sea exchange than the ocean flux by Takahashi et al. (2009) The two models showshowed better agreement for XCO2 ( Fig.23 (b) and (e)), but also differdiffered by ~1 ppmv over Taklamakan Desert and along the eastern side of the Tibet Plateau. The OCO-2 swath data arewere integrated into the corresponding horizontal grids of WRF-VPRM and CT2019, respectively, to validate XCO2. Biases of WRF-VPRM and 180 CT2019 both fall into the range of ±3 ppmv as shown in Fig.2(d3(c) and (f), respectively, but WRF-VPRM apparently providesprovided more details of spatial gradient. WRF-VPRM showsshowed well-mixed underestimations and overestimations along neighbouring satellite tracks, while CT2019 tendstended to overestimate (underestimate) over Tibet Plateau (Taklamakan Desert) where WRF-VPRM givesgave slightly smaller biases. Fig.4(a) and (b) present the raw data pairs between models and OCO-2 with daily interval for WRF-VPRM and CT2019, respectively. In general, the WRF-VPRM model 185 reproducesreproduced OCO-2 well, with mean bias (MB) of 0.76 ppmv, and normalized mean bias (NMB) of 0.19% (Fig.3(a)); CT2019 showsshowed MB of 0.54 ppmv and NMB of 0.17% (Fig.3(b)),, suggesting an overall acceptable performance of the weather-biosphere model to reproducesimulate the spatial distribution pattern of XCO2 in China.
We further analyse WRF-VPRM validation against OCO-2 for the seven vegetation types in each season and find no prominent difference (evaluation statistics summarized in Table 1). Regarding vegetation type, the model showsshowed the largest 190 normalized mean bias (NMB) MB over deciduous forest of -0.25%1.01 and 0.31%1.27 ppmv in summer and winter, respectively, both over deciduous forest which only coverscovered a very small portion in northeast China (see dominant vegetation types in Fig.1(c)).. The three most abundant coverage vegetation types in China are grass, crops, and mixed forest.
XCO2 simulated by WRF-VPRM over grass areas iswas slightly overestimated by 0.0831~0.16%68 ppmv throughout the year, and the NMBMB over mixed forest iswas -0.11%~43~0.15%,59 ppmv, indicating a good performance of the model over the 195 vast majority of areas of China. Performance over crops generally showsshowed larger discrepancy than other vegetation types, with NMBMB ranging from 0.16%66 ppmv in summer to 0.29%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-VPRM 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 200 literature reportsliteratures reported substantially different rates of ecosystem respiration and photolysis uptake (Gao et al., 2018;Yang et al., 2016;. 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 performs bestshowed slightly smaller bias in summer (NMB=0.12%) and worstlarger bias in winter (NMB=0.23%),, and the correlation coefficients arewere all ~0.8, consistent with application over the U.S. (Hu et al., 2020) which also reported slightly better performance in summer 205 than other seasons, indicating good agreement with the OCO-2 satellite product. Fig.34 also presents the overall simulation bias against ground-based observations employed in this study at thetheir raw temporal intervals (dailymonthly for OCO-2, daily for TCCONdata at HefeiESRL sites, hourly for tower data at Lin'an, and monthlydaily for dataTCCON at Hefei). At the ESRL sites). Surface (Fig.4(c)), surface CO2 concentrations arewere simulated well with minor overestimation by 0.69 ppmv (0.17%) at the ESRL sites ( Fig.3(f)). However, evaluation. Evaluation at the 210 Lin'an station shows significant overestimations for CO2 by 5.34 ppmv (1.25%) and 5.41 ppmv (1.27%) at 21m (Fig.3(d)) and 55m ( Fig.3(e)) AGL, respectively; thewas performed with the 4km-grid simulation. The mid-level heights of WRF-VPRM's first, second, and third layers arewere 12.3m, 36.9m, and 61.6m, respectively, and simulations arewere linearly interpolated to 21m and 55m to compare with the tower data. 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 215 largely attributablelikely due to the combined effect of vertical allocation of anthropogenic emission within the model as recently recognized (Brunner et al., 2019). Biosphere and parameterization of VPRM. Tracer transport models (such as WRF-VPRM and CASA) and inverse modelling methods allocate anthropogenic CO2 emission into the near surface layer due to lack of injection height information, which will likelymay subsequently lead to systematic overestimation of surface CO2 concentration in industrial areas; though. Through a regional scale (750×650km) modelling study around the city of Berlin (, 220 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-VPRM did not show significantly different performance over different vegetation types as shown in Table 1. As compared to the ESRL background sites which were located in more remote areas with little anthropogenic emission ( Fig.1(a)), Lin'an was 225 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-VPRM 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 simulations showed relatively larger bias at 55m than 21m due to smaller topography roughness and higher wind speed which 230 increases with height according to observations ( Figure S3). CT2019 also substantially overestimatesoverestimated at Lin'an, but the first, second, and third layers' mid-level heights are 25m, 103m, and 247m, respectively, so we diddo not compare it directlypresent the direct comparison with the tower data, but analysed the simulated diurnal variation as will be discussed in Section 3.3. Fig.3(c) and (d) reveal that observed average CO2 concentrations at Lin'an (428 ppmv) are substantially higher than those at ESRL sites (407-410 ppmv). The evaluation at Lin'an station also infers the prominent high CO2 level in YRD 235 due to the intensive regional anthropogenic emission as compared with ESRL sites at remote locations. Pu et al. (Pu et al., 2014)

CO2 seasonal variation and trend in China
We next analyse the seasonality of CO2 and XCO2 and find that the terrestrial flux playsplayed a more influential role than anthropogenic emission. WRF-VPRM successfully reproducesreproduced seasonal variations of CO2 at ESRL sites, with a correlation coefficient of 0.90 ( Fig.45(a)), but the )). The WRF-VPRM 4km-grid simulation showed a correlation between 250 simulated and observed CO2 at coefficient of 0.82 with the Lin'an tower is only 0.67 (Fig.4(c)); we will probe into bias at Lin'an in the next section.observation (averaged for daytime 21m and 55m data). Both the model and measurements showshowed prominent seasonal cycles for surface CO2 concentrations, with maximums. The WRF-VPRM simulation showed maxima in April (413-419 ppmv) and minimumsminima in August (399-404 ppmv) as shownpresented in Fig.45(b). The model suggestssuggested that the anthropogenic CO2 contribution iswas 2.6 ppmv in both months, while the biogenic 255 contributions arewere 3.1 ppmv and -1.2 ppmv in April and August, respectively (Fig.45(d)). Anthropogenic emission ( Fig.45(f)) showsshowed a flat curve with relatively higher values in December due to fuel combustion for heating (Zheng et al., 2018);. EVI meanwhile shows maximumsshowed maxima in July and August (Fig.45(f)). During summer, photosynthetic uptake almost completely offsetscompensated anthropogenic emission, causing the minimum CO2 concentration observed in August, while the higher anthropogenic emission in December and respiration flux in April leadled to the two corresponding 260 peaks. The anthropogenic XCO2 contributions from arewere 0.5 and 0.6 ppmv in April and August, respectively, and the biogenic contributions arewere 0.8 ppmv and -1.5 ppmv, respectively, suggesting that the seasonality of XCO2 iswas also primarily dominated by terrestrial flux. Furthermore, the seasonality at high-latitude ESRL sites (UUM and WLG) iswas stronger than at Lin'an and low-latitude sites (DSI and LLN) because of the larger temperature and photosynthetically active radiation (PAR) gradients. Annual average anthropogenic and biogenic XCO2 contributions arewere 7.1 ppmv and -1.9 ppmv, 265 respectively, indicating that biosphere uptake iswas an important carbon sink offsetting 27% of anthropogenic emission and slowing the growth of atmospheric CO2.
XCO2 showsshowed similar seasonality, with minimumsminima in August and maximumsmaxima in April and December ( Fig.45(b)). Both WRF-VPRM and CT2019 showshowed good agreement with TCCON Hefei observations with correlations of 0.89 and 0.88, respectively (Fig.45(e)). However, we note that WRF-VPRM simulatessimulated drastic changes (e.g., the 270 grey shaded period in Fig.45(e)) that arewere not shownreproduced by CT2019;. Fig. 56 shows the daily concentrations of XCO2 overlaid with horizontal wind speed at 10m AGL from WRF-VPRM and CT2019 and highlights large discrepancies over Hefei. ( Figure S4 shows the same comparison but using WRF-VPRM 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 shown in Fig.56(g)-(i)). At the leading edge of the front, a convergence zone associated with a low 275 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 hashad difficulty reproducing such regional weather systems that can lead to rapid and localized changes in CO2 concentration and terrestrial 280 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 is 0.56%/yr (was 2.2 ppmv/yr (0.56%/yr) at the ESRL sites and 2.3 ppmv/yr (0.67%/yr (2.8 ppmv/54%/yr) at Lin'an. The observed average CO2 concentrations at Lin'an (428 ppmv) were substantially higher than those at ESRL sites (407-410 ppmv). The prominent higher 285 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 iswas also found to increase by 0.57%/yr (2.3 ppmv/yr (0.57%/yr) as suggested by OCO-2 and 0.61%/yr (2.5 ppmv/yr (0.61%/yr) as suggested by the simulation. WRF-VPRM reproducesreproduced the trends in good agreement with ground and satellite observations.
Model simulated budgets suggestssuggested that the increasing trends for anthropogenic, biogenic, and background XCO2 290 arewere 0.81%/yr, -9.17%/yr, and 0.59%/yr, respectively; the trends for anthropogenic, biogenic, and background CO2 arewere 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)  in East Asia, photosynthetic uptake also servedserves 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 China. 300

Diurnal variation of near-surface CO2 and influence factors
Finally, we examine the diurnal variation of meteorology and CO2 data at Lin'an station as shown in Fig.67 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 arewere larger in summer (JJA) than winter (DJF) and arewere larger at 21m than 55m, indicating the dominant influence of terrestrial flux over anthropogenic emission in determining the near surface CO2 305 concentration. Fig.67(c) presents the diurnal change of wind speed collected at 50m of the Lin'an tower. The higher wind speed between 10:00-22:00 local time suggests strong regional transport and mixing of CO2 mainly occurs during this period. Fig.6(d) and (g) present the WRF-VPRM and CT2019 simulation bias, at 21m and 55m respectively, against Lin'an tower data at 21m (note the Y-scales are different).and Fig.7(e) and (f) present the bias of CT2019 at 21m and 55m respectively. We find that both models prominently overestimate during night time.overestimated during nighttime, which shall be attributed to 310 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 ECeddy-covariance tower measurement.
Fig.6(f7(g) and (ih) present the simulated NEE by WRF-VPRM 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-VPRM may have also estimatedoverestimate nighttime ecosystem respiration during the non-growing season,at Lin'an as it has a warmer climate 315 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 atmospheric CO2 as shown in Fig.78(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, 320 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., 2020;Hu et al., 2012).
Such boundary layer characteristics are confirmed by CO2 vertical gradients at Lin'an. in this study. CO2 at 55m height iswas 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 depletesdepleted the near surface CO2 and daytime boundary layer convection 325 developsdeveloped, the CO2 gradient iswas gradually weakened from 06:00 to 11:00 LT and remainsremained minimal through the rest of the daytime; at midday when photosynthesis reaches maximum intensity, CO2 at 21m iswas even lower than at 55m. WRF-VPRM generally reproducesroughly reproduced the diurnal profile but noticeably underestimatesunderestimated the intensity of night time 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). 330 The relationship between the near-surface CO2 profile and boundary layer stability is further statistically examined. Fig.78(b) presents the correlation between air temperature gradient (ΔT/ΔH) and CO2 concentration gradient (ΔCO2/ΔH) calculated with annual averaged diurnal profiles of tower observations, which 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 335 Δ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 Fig.78(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.78(b)) implying a lower CO2 concentration at the surface. While the diurnal variations of ΔCO2 340 arewere 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.78(a)) likely suggest influence of transport of CO2 from urban plumes in the region; for example, from Hangzhou which is 60 km away from the tower. Due to rush-hours anthropogenic emissions, CO2 enhancement at Hangzhou relative to a background site exhibitsexhibited 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 345 cities such as Hangzhou may be transported to the Lin'an site. Due to higher altitude and stronger windswind profile increases with height at Lin'an according to observations ( Figure S4) -55m at the Lin'an tower has a larger footprint than 21m, thus 55m on the tower is more likely affected by the urban plumes in the region than 21m. The 10:00 and 18:00 LT ΔCO2 peaks at Lin'an likely suggest stronger CO2 enhancement at 55m than at 21m from influence of regional anthropogenic emissions; the slight delay of these ΔCO2 peaks relative to rush hours (at about 08:00 and 17:00 LT) further corroborate the 350 hypothesis of delayed influence of transport of urban CO2 from Hangzhou. Even though 55m has a larger footprint than at 21m and thus may be more likely affected by regional urban emissions, turbulent vertical mixing may reduce the different influence from regional urban emissions, which explains the fact that ΔCO2 peaks are only minor. 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), but our study appliesapplied tower data as direct evidence 355 to demonstrate the significant impact of PBL thermal structure, which has rarely been documented. More importantly, although WRF-VPRM failsfailed to capture the bimodal ΔCO2 peaks at rush hours, because monthly ODIAC data lackslacked 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.

Summary and Conclusions 360
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 "background" stations, indicating the dominant influence of anthropogenic emission at regional scales. The Lin'an tower data 365 shows a large negative correlation (-0.8598 A steadily increasing trend in XCO2 by ~2.46 ppmv (~0.6%/yr) during the study period is demonstrated consistently by both 375 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 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 380 the world.
The most significant modelling bias is identified from validation against the Lin'an tower data55m observations, which WRF-VPRM 4km-grid simulation overestimated by about 5.381.06 ppmv (1.26%) with a correlation coefficient of 0.6782. The allocation of anthropogenic emission into the surface layer is partially responsible for this modelling bias because Lin'an is closely affected by upwind industrial mega cities in YRD, suggesting the need to include vertical profiles of fossil fuel 385 combustion to properly redistribute the ODIAC for modelling purposes. HoweverIn addition, diurnal variations of the bias suggest that the modelling discrepancy is likely due toalso induced by large uncertainty associated with simulating nighttime ecosystem respiration during the nighttime. Representation and parameterization of photosynthetic carbon uptake in VPRM has been continuously improved during the past 10 years since its first release (Hu et al., 2020), but ecosystem respiration parameterization is still too simplified to fully represent the autotrophic and heterotrophic respiration of biomass. Li et al. 390 (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-VPRM 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 simultaneously, it promotes the opportunity to investigate the interactions between atmospheric mixing 395 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 shortterm carbon cycles at regional scales.

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

Author contributions
The concept and ideas to design the integrated simulation and observation analysis wereare devised by YJ, X-MH, and XD. 405 Simulation wasis performed by X-MH. OCO-2 satellite product wasis collected and processed by X-MH. CT2019 assimilation data and ground-based observations wereare collected by XD. Tower measurement wasis conducted, processed, and analysed by QM, JP, and YJ. Model evaluation wasis performed by MY. The manuscript wasis prepared by XD and X-MH with input and feedback from YJ, MY, QM, JP, and GZ. 410 Competing interests The authors declare that they have no conflict of interest.

Acknowledgements
This work is supported by the Fundamental Research Funds for the Central Universities (14380049)

650
The blue box represents location of Hefei.