Analysis of CO₂ spatio-temporal variations in China using a weather–biosphere online coupled model

The dynamics of atmospheric CO₂ has received considerable attention in the literature, yet significant uncertainties remain within the estimates of contribution from the terrestrial flux and the influence of atmospheric mixing. In this study we apply the WRF-Chem model configured with the Vegetation Photosynthesis and Respiration Model (VPRM) option for biomass fluxes in China to characterize the dynamics of CO₂ in the atmosphere. The online coupled WRF-Chem model is able to simulate biosphere processes (photosynthetic uptake and ecosystem respiration) and meteorology in one coordinate system. We apply WRF-Chem for a multi-year simulation (2016–2018) with integrated data from a satellite product, flask samplings, and tower measurements to diagnose the spatio-temporal variations of CO₂ fluxes and concentrations in China. We find that the spatial distribution of CO₂ was dominated by anthropogenic emissions, while its seasonality (with maxima in April 15 ppmv higher than minima in August) was dominated by the terrestrial flux and background CO₂. Observations and simulations revealed a consistent increasing trend in column-averaged CO₂ (XCO₂) of 2.46 ppmv (0.6 % yr⁻¹) resulting from anthropogenic emission growth and biosphere uptake. WRF-Chem successfully reproduced ground-based measurements of surface CO₂ concentration with a mean bias of −0.79 ppmv and satellite-derived XCO₂ with a 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 at Lin'an (30.30^∘ N, 119.75^∘ E) revealed a strong correlation (−0.98) between vertical CO₂ and temperature gradients, suggesting a significant influence of boundary layer thermal structure on the accumulation and depletion of atmospheric CO₂.

1 Introduction

Climate research requires accurate characterization of atmospheric CO₂, which is closely affected by both atmospheric transport and terrestrial sources and sinks (Bauska et al., 2015; Keenan et al., 2016). Our current knowledge largely comes from interpreting ground- or space-based measurements and model simulations. While observation is limited by spatial and temporal coverage, modelling approaches also suffer from various uncertainties (Shi et al., 2018). Modelling assessment of CO₂ 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); second, inverse modelling in which prior flux estimates applied in atmospheric transport models are adjusted by observational data and/or satellite products to determine the 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 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 estimates of prior flux (Crowell et al., 2018; Feng et al., 2019), and advanced measurement and 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 the 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. Previous modelling studies (Ahmadov et al., 2009; Kretschmer et al., 2012; Park et al., 2018, 2020; Beck et al., 2013; Pillai et al., 2012) have demonstrated the weather–biosphere coupled model can successfully capture the mesoscale CO₂ transport at regional and local scales with significant improvements. But whether it can reproduce the long-term variations and subsequently estimate carbon fluxes at regional scales with high confidence remains a crucial issue to be addressed.

Understanding the spatio-temporal characteristics of atmospheric CO₂ 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 CO₂ at regional to national scales remains poorly understood due to a 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 CO₂ fluxes and concentrations over mixed forest and rice paddy areas in northeast China, but the 1-year simulation limits the attempt to investigate interannual CO₂ variation, which is subject to substantial change (Z. Fu et al., 2019). 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 the global terrestrial flux. Fu et al. (2019) applied the GEOS-Chem simulation with the offline Carbon Tracker (Peters et al., 2007) as input to estimate the impacts of the terrestrial flux and anthropogenic emissions on the annual variation of CO₂ concentrations, but regional-scale assessment was limited by the coarse grid resolution (2^∘ × 2.5^∘). A machine learning technique has also been employed to upscale site observations to the 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 the coarse grid resolution. These pilot studies have shed light on improving the understanding of spatio-temporal characteristics of CO₂ in China with modelling or observational methods, but an integrated investigation with both modelling and observations at a fine scale is urgently needed to expand diagnostic understanding of localized and regional transport, flux, and concentration of CO₂ to inform emission management and climate adaption policies (Y. Fu et al., 2019; Niu et al., 2017; Wang et al., 2019).

In this study we use the WRF-Chem model configured with the Vegetation Photosynthesis and Respiration Model (VPRM) option (Hu et al., 2020; Mahadevan et al., 2008) to simulate and characterize the spatio-temporal variation of atmospheric CO₂ in China from 2016–2018 and also to validate this weather–biosphere model with recent advanced satellite and tower observations. WRF-Chem 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 attempts to apply and evaluate it for a multi-year simulation over China. We first describe the modelling methods and data followed by model validation against observations from multiple datasets and then present the spatio-temporal variations and estimates of contributions from anthropogenic emissions, the terrestrial flux, and background concentrations. Finally, we investigate tower data and reveal the boundary layer thermal structure impacts on atmospheric CO₂ accumulation and depletion.

Figure 1 Annual averages of (a) ODIAC emission, (b) MODIS EVI, and (c) dominant vegetation type in the 20 km simulation domain and (d) terrain height of the 4 km simulation domain. The locations of the ESRL sites, TCCON Hefei site, and Lin'an tower site are indicated with red circles, rectangles, and diamonds, respectively in (a). The 4 km domain is indicated with the dashed red rectangle in (b), and the locations of Hangzhou and Shanghai are indicated with yellow triangles in (d).

2 Method

We conduct nested WRF-Chem (Version 3.9.1.1) simulations over China (domain shown in Fig. 1a) and the Yangtze River Delta (YRD) region (domain shown in Fig. 1d) at 20 and 4 km grid resolution, respectively. Both simulations were configured with 47 vertical layers with model tops at 10 hPa. 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., 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 the Supplement. In general, the 4 km grid simulation showed no significant difference as compared to the 20 km grid simulation (demonstrated in Figs. S1 and S2 in the Supplement); thus the 20 km grid simulation was used to characterize the spatio-temporal distributions of CO₂ over China, and the 4 km grid simulation was only used for comparison with tower data collected at a background site in YRD. Discussions in the next section will mostly refer to the 20 km grid simulation unless otherwise specified. Initial and lateral boundary conditions for the 20 km grid simulations were derived from the mole fraction product of CarbonTracker (Peters et al., 2007) with 3^∘ × 2^∘ resolution. The latest update of column-averaged CO₂ (XCO₂) concentration assimilation product from CarbonTracker (CT2019) with 1^∘ × 1^∘ resolution (Jacobson, 2020) was also employed for comparison with the WRF-Chem simulation. The anthropogenic emission inventory is the Open-Data Inventory for Anthropogenic Carbon dioxide (ODIAC) with 0.1^∘ × 0.1^∘ resolution (Oda et al., 2018), shown in Fig. 1a. ODIAC has been widely applied in recent modelling studies and demonstrated good agreement with other global inventories (Hedelius et al., 2017; Hu et al., 2020). The ocean flux is from climatology estimation (Takahashi et al., 2009); and vegetation fractions and enhanced vegetation index (EVI; shown in Fig. 1b) are from MODIS (Huete et al., 2002). CO₂ from initial and boundary conditions, anthropogenic emission, and the terrestrial biogenic flux was tagged as BCG, ANT, and BIO, respectively, to allow the contributions from each process to be identified and quantified through one simulation.

WRF-Chem calculates ecosystem respiration (ER) and gross ecosystem exchange (GEE) with the following functions as

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{ER}=\mathit{\alpha }\cdot T+\mathit{\beta }\text{(2)}& {\displaystyle }\begin{array}{rl}\text{GEE}=& -\mathit{\lambda }\cdot T_{\text{scale}}\cdot W_{\text{scale}}\cdot P_{\text{scale}}\cdot {\left(\mathrm{1}+\text{PAR}/{\text{PAR}}_{\mathrm{0}}\right)}^{-\mathrm{1}}\\ & \cdot \text{EVI}\cdot \text{PAR},\end{array}\end{array}$$

where T is the air temperature at 2 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (T2); α, β, and λ are vegetation-type-dependent parameters; PAR₀ is the vegetation-type-dependent half-saturation value of photosynthetically active radiation (PAR); and T_scale, W_scale, and P_scale 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-Chem is parameterized (α, β, λ) for seven vegetation types (Fig. 1c): 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 eddy-covariance (EC) tower measurements 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 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 (R_H). NPP is the difference between photosynthetic uptake (equivalent to GEE) and autotrophic respiration (R_A). The summary of R_H and R_A is equivalent to ER. Thus, WRF-Chem 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 R_A, 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 more details in Sect. 3.2.

Figure 2 Photos of the (a) Lin'an regional atmospheric background station and (b) the data analysis lab. Wind rose map at Lin'an derived from wind speed and wind direction observations for 2016–2018 at (c) 10 m and (d) 50 m.

Hourly measurements of CO₂ concentrations were collected at the Lin'an Regional Atmospheric Background Station (30.30^∘ N, 119.75^∘ E, surroundings shown in Fig. 2a) with Picarro G1301 and G1302 trace gas analysers mounted on an observation tower at 21 and 55 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ (above ground level), respectively, and analysed online (data analysis lab shown in Fig. 2b). The station is located in the remote area of Hangzhou 138.6 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in the middle of a hilly area covered by mixed forest. The observation tower is 60 km to the west of the downtown centre of Hangzhou and 195 km to the southwest of Shanghai. Figure 2c and d present the wind rose map at Lin'an derived from hourly observations of 10 and 55 m 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 (Deng et al., 2018; Pu et al., 2020). The tower data provide a representative sampling of CO₂ gradients resulting from the exchange between atmosphere mixing and the terrestrial flux.

Atmospheric samples near the surface were collected at monthly intervals and analysed for CO₂ through the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) at four sites (locations shown in Fig. 1a) within our study domain, including Dongsha Island (DSI; 20.69^∘ N, 116.73^∘ E), Lulin (LLN; 23.47^∘ N, 120.87^∘ E), Ulaan Uul (UUM; 44.45^∘ N, 111.09^∘ E), and Mt. Waliguan (WLG; 36.29^∘ N, 100.89^∘ E). The Orbiting Carbon Observatory-2 (OCO-2) satellite product (Kiel et al., 2019) with daily intervals was employed to validate simulation of column-averaged CO₂ (XCO₂) concentrations. A total of 204 940 OCO-2 version 9 swath data covering the simulation period were used in this study. Daily ground-based Fourier transform spectrometer (FTS)-measured XCO₂ data at Hefei site (31.90^∘ N, 117.17^∘ E) were also collected through the Total Carbon Column Observing Network (TCCON) for the year 2016 (Wang et al., 2017). The TCCON Hefei site was located in the northwestern rural area of the city of Hefei, and measurements were conducted from September 2015 to December 2016 (Wang et al., 2017). 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 CO₂, which has not been thoroughly discussed in the literature.

Figure 3 2016–2018 averages of WRF-Chem simulations of (a) first-layer (mid-level height is 12 km) CO₂ concentration and (b) XCO₂ concentration. (c) WRF-Chem-simulated XCO₂ bias against OCO-2. Panels (d)–(f) are the same as (a)–(c) but for CT2019 (first-layer mid-level height is 25 m).

Figure 4 Data pairs for OCO-2 against (a) WRF-Chem and (b) CT2019; (c) ESRL against WRF-Chem; Lin'an tower against WRF-Chem 4 km grid simulation at (d) 21 m and (e) 55 m; and (f) TCCON Hefei against WRF-Chem and CT2019.

3 Result and discussion

3.1 Model evaluation

We first evaluate the capability of WRF-Chem to reproduce concentrations of surface CO₂ and XCO₂, and we find fairly good model performance through the comparison with satellite and ground-based observations. The WRF-Chem-simulated surface layer (mid-level height a.g.l. is 12 m) CO₂ and XCO₂ averages between 2016–2018 are demonstrated in Fig. 3a and b, respectively. High concentrations were found over industrial areas such as the North China Plain (NCP), Pearl River Delta (PRD), and Yangtze River Delta (YRD), where the surface CO₂ and XCO₂ were above 440 and 408 ppmv, respectively; the domain averages were 411 and 406 ppmv, respectively. While most climate models assume evenly distributed CO₂ (Fung et al., 1983; Kiehl and Ramanathan, 1983), our data demonstrate a prominent gradient between industrial and remote areas (e.g., Tibetan Plateau, Mongolia), especially for surface CO₂, which could be an important source of uncertainty for estimating the long-wave radiation effect (Xie et al., 2018). Spatial patterns of CO₂ and XCO₂ were in close agreement with ODIAC, indicating the dominant impact of anthropogenic emission in determining the CO₂ distribution. WRF-Chem-simulated CO₂ was generally consistent with CT2019 (Fig. 3c), but CT2019 estimated near-surface CO₂ (mid-level height a.g.l. is 25 m) 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-Chem. The two models showed better agreement for XCO₂ (Fig. 3b and e) but also differed by ∼ 1 ppmv over Taklamakan Desert and along the eastern side of the Tibetan Plateau. The OCO-2 swath data were integrated into the corresponding horizontal grids of WRF-Chem and CT2019, respectively, to validate XCO₂. Biases of WRF-Chem and CT2019 both fall into the range of ± 3 ppmv as shown in Fig. 3c and f, respectively, but WRF-Chem evidently provided more details of the spatial gradient. WRF-Chem showed well-mixed underestimations and overestimations along neighbouring satellite tracks, while CT2019 tended to overestimate (underestimate) over the Tibetan Plateau (Taklamakan Desert) where WRF-Chem gave slightly smaller biases. Figure 4a and b present the raw data pairs between models and OCO-2 with daily intervals for WRF-Chem and CT2019, respectively. In general, the WRF-Chem model reproduced OCO-2 well, with a mean bias (MB) of 0.76 ppmv, and CT2019 showed a MB of 0.54 ppmv, suggesting an overall acceptable performance of the weather–biosphere model to simulate the spatial distribution pattern of XCO₂ in China.

We further analyse WRF-Chem 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 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. XCO₂ simulated by WRF-Chem 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 China. Performance over crops generally showed a larger discrepancy than other vegetation types, with the MB ranging from 0.66 ppmv in summer to 1.19 ppmv in winter, suggesting the model tends to slightly overestimate column concentration of CO₂ over cropland. Li et al. (2020) reported that WRF-Chem underestimated biosphere carbon over rice paddy sites (by ∼ 3 %) in northeast China and suggested the parameterization of α, β, and λ as the most important cause. Cropland differs significantly across China, with various types of species such as rice, wheat, and corn, for which literature reported 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 XCO₂. In terms of seasonal difference, WRF-VPRM showed a slightly smaller bias in summer and a larger bias in winter, and the correlation coefficients were all ∼ 0.8, consistent with application over the United States (Hu et al., 2020), which also reported slightly better performance in summer than other seasons.

Figure 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. 4c), surface CO₂ concentrations were simulated well with minor overestimation by 0.69 ppmv. Evaluation at the Lin'an station was performed with the 4 km grid simulation. The mid-level heights of WRF-Chem's first, second, and third layers were 12.3, 36.9, and 61.6 m, respectively, and simulations were linearly interpolated to 21 and 55 m for comparison with the tower data. The evaluation at 21 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ (Fig. 4d) shows slight overestimation by 0.02 ppmv, but the evaluation at 55 m height (Fig. 4e) 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-Chem and CASA) and inverse modelling methods allocate anthropogenic CO₂ emission to the near-surface layer due to lack of injection height information, which may subsequently lead to systematic overestimation of surface CO₂ concentration in industrial areas. Through a regional-scale (750 km × 650 km) modelling study around the city of Berlin, Brunner et al. (2019) reported that distributing anthropogenic emission into the surface layer caused near-surface CO₂ concentration to be overestimated 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 in a densely vegetated area. Validation against OCO-2 suggested that WRF-Chem 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. 1a), 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 a more prominent error at Lin'an. In fact, the 20 km grid WRF-Chem simulation bias at Lin'an was 5.34 and 5.41 ppmv at 21 and at 55 m, respectively (Fig. S2), significantly larger than the bias at ESRL sites. In addition, both the 20 km grid and 4 km grid simulations showed relatively larger bias at 55 m than 21 m due to smaller topography roughness and higher wind speed, which increases with height according to observations (Fig. S3). CT2019 also substantially overestimated at Lin'an, but the first, second, and third layers' mid-level heights are 25, 103, and 247 m, respectively, so we do not present a direct comparison with the tower data. Simulated XCO₂ from both WRF-Chem and CT2019 was consistent with the TCCON Hefei site observations as shown in Fig. 4f, with a MB of −0.79 and −0.78 ppmv, respectively, and a normalized mean bias (NMB) of −0.20 % and −0.19 %, respectively. The 4 km grid simulation showed a very similar result to the 20 km grid simulation for XCO₂ (Figs. S1 and S2). Recent atmospheric inverse modelling studies (Y. Fu et al., 2019; Wang et al., 2019; Xie et al., 2018) reported the simulation bias of XCO₂ as 0.5–2 ppmv with posterior flux inputs. The WRF-Chem model applied in this study has demonstrated good agreement with the observations though our evaluation.

Figure 5 Monthly variations of (a) CO₂ at ESRL sites, (b) total (black) and background (BCG, grey) CO₂ (line) and XCO₂ (area and bar), (c) CO₂ at Lin'an station (averaged for daytime 21 and 55 m data). (d) Contributions from anthropogenic (ANT, orange) and biogenic (BIO, blue) for CO₂ (lines) and XCO₂ (bars); (f) ODIAC emission and MODIS EVI; and (e) daily variation of XCO₂ at the TCCON Hefei site.

3.2 CO₂ seasonal variation and trend in China

We next analyse the seasonality of CO₂ and XCO₂ and find that the terrestrial flux played a more influential role than anthropogenic emission. WRF-Chem successfully reproduced seasonal variations of CO₂ at ESRL sites, with a correlation coefficient of 0.90 (Fig. 5a). The WRF-Chem 4 km grid simulation showed a correlation coefficient of 0.82 with the Lin'an tower observation (averaged for daytime 21 and 55 m data). Both the model and measurements showed prominent seasonal cycles for surface CO₂ concentrations. The WRF-Chem simulation showed maxima in April (413–419 ppmv) and minima in August (399–404 ppmv) as presented in Fig. 5b. The model suggested that the anthropogenic CO₂ 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. 5d). Anthropogenic emission (Fig. 5f) 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. 5f). During summer, photosynthetic uptake almost completely compensated for anthropogenic emission, causing the minimum CO₂ concentration to be observed in August, while the higher anthropogenic emission in December and respiration flux in April led to the two corresponding peaks. The anthropogenic XCO₂ 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 XCO₂ was also primarily dominated by the 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) gradients. Annual average anthropogenic and biogenic XCO₂ 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 CO₂.

Figure 6 Daily XCO₂ from CT2019 (a–c) and WRF-Chem (d–f). Weather map from the Korea Meteorological Administration (g–i). The blue box represents the location of Hefei.

Table 1 Evaluation statistics^a for WRF-Chem 20 km grid simulation against the OCO-2 satellite product at daily intervals.

^a Mean bias was calculated as $\text{MB}=\frac{\mathrm{1}}{N}\sum _{i=\mathrm{1}}^N({\text{Sim}}_i-{\text{Obs}}_i)$, and correlation coefficient was calculated as $\text{cc}=\frac{\sum _{i=\mathrm{1}}^N\left({\text{Sim}}_i-\stackrel{\mathrm{‾}}{\text{Sim}}\right)\left({\text{Obs}}_i-\stackrel{\mathrm{‾}}{\text{Obs}}\right)}{\sqrt{\sum _{i=\mathrm{1}}^N{\left({\text{Sim}}_i-\stackrel{\mathrm{‾}}{\text{Sim}}\right)}^{\mathrm{2}}}\sqrt{\sum _{i=\mathrm{1}}^N({\text{Obs}}_i-\stackrel{\mathrm{‾}}{\text{Obs}}{)}^{\mathrm{2}}}}$, where $\stackrel{\mathrm{‾}}{\text{Sim}}$ is the average of simulations, and $\stackrel{\mathrm{‾}}{\text{Obs}}$ is the average of observations. ^b For each season, the evaluation statistic with the worst performance (largest absolute value of MB) is highlighted in italics, and the one with the best performance (smallest absolute value of MB) is highlighted in bold.

XCO₂ showed similar seasonality, with minima in August and maxima in April and December (Fig. 5b). Both WRF-Chem and CT2019 showed good agreement with TCCON Hefei observations, with correlations of 0.89 and 0.88, respectively (Fig. 5e). However, we note that WRF-Chem simulated drastic changes (e.g., the grey shaded period in Fig. 5e) that were not reproduced by CT2019. Figure 6 shows the daily concentrations of XCO₂ overlaid with horizontal wind speed at 10 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ from WRF-Chem and CT2019 and highlights large discrepancies over Hefei (Fig. S4 shows the same comparison but using WRF-Chem 4 km grid simulation data). Between 1 and 3 April 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. 6g–i). At the leading edge of the front, a convergence zone associated with a low-pressure centre 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 CO₂. Meanwhile, the cold front transported anthropogenic CO₂ from NCP to YRD, and the convergence zone along YRD ahead of the front facilitated the accumulation of air pollutants and CO₂ from anthropogenic emissions. With its coarse spatio-temporal resolution, CT2019 had difficulty reproducing such regional weather systems that can lead to rapid and localized changes in CO₂ concentration and the terrestrial flux, indicating the importance of fine-resolution modelling to better represent the small spatial-scale and rapid temporal-scale variations of CO₂ (Agusti-Panareda et al., 2019).

We also find a notable increasing trend for the 3-year study period. Observed CO₂ 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 CO₂ concentrations at Lin'an (428 ppmv) were substantially higher than those at ESRL sites (407–410 ppmv). The prominent higher levels of CO₂ 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 a rate of 0.82 % yr⁻¹ as suggested by ODIAC. Domain-wide XCO₂ 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-Chem reproduced the trends in good agreement with ground and satellite observations. Model-simulated budgets suggested that the increasing trends for anthropogenic, biogenic, and background XCO₂ were 0.81, −9.17, and 0.59 % yr⁻¹, respectively; the trends for anthropogenic, biogenic, and background CO₂ were 4.95, −0.73, 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 the terrestrial flux. Wu et al. (2012) reported that measured CO₂ concentration at the Changbai Mountain forest site in northeast China increased by 1.76 ppmv yr⁻¹ between 2003 and 2010. With the atmospheric inversion modelling method, Z. Fu et al. (2019) estimated that surface CO₂ in East Asia increased by 2–3 ppmv yr⁻¹ between 2004 and 2012. These trends suggest that 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 the terrestrial flux is usually critically large and is affected by both vegetation itself and climate conditions (Z. Fu et al., 2019; Niu et al., 2017), simulation over longer time periods is necessary in future studies to conclusively comment on the changing trend of biosphere CO₂ in China.

Figure 7 Seasonal mean diurnal variations of observed CO₂ at (a) 21 m and (b) 55 m. WRF-Chem simulation biases of CO₂ at (c) 21 m and (d) 55 m. CT2019-simulated biases at (e) 21 m and (f) 55 m. Simulated NEE from (g) WRF-Chem and (h) CT2019.

3.3 Diurnal variation of near-surface CO₂ and influence factors

Finally, we examine the diurnal variation of CO₂ data at Lin'an station as shown in Fig. 7 to reveal the temporal dynamics and atmospheric mixing of CO₂ at a local scale. While both 21 m (Fig. 7a) and 55 m (Fig. 7b) CO₂ show prominent diurnal changes, the variations were larger in summer (JJA) than winter (DJF) and were larger at 21 m than 55 m, indicating the dominant influence of the terrestrial flux over anthropogenic emission in determining the near-surface CO₂ concentration. Figure 7c and d present the WRF-Chem simulation bias at 21 and 55 m, respectively, and Fig. 7e and f present the bias of CT2019 at 21 and 55 m, respectively. We find that both models prominently overestimated during night-time, which shall be attributed to the bias in simulating NEE. Li et al. (2020) reported that the model overestimated night-time NEE at a mixed forest site, Wuying (47.15^∘ N, 131.94^∘ E), by 34 % during the growing season (May–September) according to eddy-covariance tower measurements. Figure 7g and h present the simulated NEE by WRF-Chem and CT2019, respectively, which show close correlations with the CO₂ simulation biases. While Lin'an is also covered by mixed forest, our evaluation suggests that WRF-Chem may also overestimate night-time ecosystem respiration at Lin'an as it has warmer climate conditions than Wuying (Fig. S5), and CT2019 has an even greater bias for presenting the diurnal cycles of the terrestrial flux.

We also find that planetary boundary layer height (PBLH) significantly affects diurnal accumulation and depletion of atmospheric CO₂ as shown in Fig. 8a. During daytime in the growing season, photosynthetic uptake results in lower CO₂ concentration; meanwhile, PBLH is also high and allows rapid vertical mixing between near-surface and upper air. During night-time when photosynthesis stops, CO₂ 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., 2020, 2012). Such boundary layer characteristics are confirmed by CO₂ vertical gradients at Lin'an in this study. CO₂ at 55 m height was consistently lower than the near-surface air at 21 m during night-time due to accumulation of respired CO₂ in the stable boundary layer. As photosynthetic uptake depleted the near-surface CO₂, and daytime boundary layer convection developed, the CO₂ gradient was gradually weakened from 06:00 to 11:00 LT and remained minimal through the rest of the day; at midday when photosynthesis reaches maximum intensity, CO₂ at 21 m was even lower than at 55 m. WRF-Chem roughly reproduced the diurnal profile but noticeably underestimated the intensity of night-time CO₂ difference (ΔCO₂), likely due to the bias for simulating the night-time terrestrial flux as discussed above or the underestimation of night-time boundary layer stability by the PBL scheme (Hu et al., 2012).

Figure 8 (a) Average (2016–2018) diurnal variations of simulated (black line) and observed (red line) ΔCO₂ and simulated (blue line) PBLH at Lin'an station and (b) correlation between CO₂ gradient between 55 and 21 m ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/\mathrm{\Delta }H$) and temperature gradient ($\mathrm{\Delta }T/\mathrm{\Delta }H$) at Lin'an station.

The relationship between the near-surface CO₂ profile and boundary layer stability is statistically examined further. Figure 8b presents the correlation between air temperature gradient ($\mathrm{\Delta }T/\mathrm{\Delta }H$) and CO₂ concentration gradient ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/\mathrm{\Delta }H$) calculated with diurnal profiles of tower observations averaged for 2016–2018, where ΔT, ΔH, and ΔCO₂ are the differences of temperature, height, and CO₂ concentration between the two tower levels, respectively. Figure 8b clearly demonstrates the influence of boundary layer stability on the CO₂ vertical profile, as the correlation between $\mathrm{\Delta }T/\mathrm{\Delta }H$ and $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/\mathrm{\Delta }H$ reaches −0.98. On the one hand, a more stable PBL with a strongly positive temperature gradient would promote surface CO₂ accumulation and lead to a strongly negative CO₂ gradient, especially under inversion conditions when upper air has higher temperature (orange area in Fig. 8b). Conversely, a strongly negative temperature gradient indicates stronger radiation and subsequently greater photosynthesis and CO₂ depletion in the near-surface layer, which would result in a positive CO₂ gradient (green area in Fig. 8b), implying a lower CO₂ concentration at the surface. While the diurnal variations of ΔCO₂ were primarily dictated by local biogenic CO₂ fluxes and boundary layer dynamics, the two minor daytime peaks of ΔCO₂ at Lin'an, at 10:00 and 18:00 LT (Fig. 8a), likely suggest influence of transport of CO₂ from urban plumes in the region, for example, from Hangzhou which is 60 km away from the tower. Due to rush-hour anthropogenic emissions, CO₂ 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 CO₂ plumes from cities such as Hangzhou may be transported to the Lin'an site. The influence of boundary layer conditions on CO₂ variability has been discussed in several studies through the analysis of mountain site ground-based observations (Arrillaga et al., 2019; Esteki et al., 2017; Li et al., 2014), 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-Chem failed to capture the bimodal ΔCO₂ peaks at rush hour because monthly ODIAC data lacked an hourly profile, our analysis reveals the critical importance of careful configuration of the PBL scheme and spatio-temporal distribution of anthropogenic emission for weather–biosphere modelling of atmospheric CO₂.

4 Summary and conclusions

In this study, the spatio-temporal variations of CO₂ 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 CO₂ 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 show a large negative correlation (−0.98) between vertical CO₂ concentration and air temperature gradients, suggesting the significant influence of boundary layer stability on CO₂ accumulation and depletion. The online coupled weather–biosphere model WRF-Chem enables process-based estimations of contributions from anthropogenic emission (0.59 ppmv, 0.15 %), the terrestrial flux (0.16 ppmv, −0.04 %), and background concentration (405.70 ppmv, 99.89 %) to average total XCO₂. Respective simulation biases of surface CO₂ and XCO₂ are 0.69 and 0.76 ppmv against ESRL site observations and the OCO-2 satellite product, with correlations of 0.87 and 0.90, indicating overall good performance of the WRF-Chem model. Maximum CO₂ concentrations are found in April, and minima are found in August for all 3 years, and the seasonality is reproduced well by the model, which also reveals that the terrestrial flux and background concentration dominated the seasonality rather than anthropogenic emission.

A steadily increasing trend in XCO₂ by 2.46 ppmv (∼ 0.6 % yr⁻¹) during the study period is demonstrated consistently by both model simulation and the satellite product. Budget analysis suggests that anthropogenic emission increased by 0.83 % yr⁻¹, contributing to the 0.81 % yr⁻¹ growth rate of anthropogenic XCO₂ enhancement, 27 % of which was offset by biosphere uptake. It is noted that the 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 XCO₂, representing contributions from global circulation, increased by 2.37 ppm (0.59 % yr⁻¹), suggesting that the CO₂ 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 55 m observations, which WRF-Chem 4 km 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 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 night-time ecosystem respiration. Representation and parameterization of photosynthetic carbon uptake in VPRM have 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 (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 are necessary to further validate the anthropogenic fluxes from ODIAC. Nevertheless, WRF-Chem is demonstrated to be a reliable tool to model the dynamics of CO₂ and exchange between the atmosphere and the 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 and the terrestrial flux (Carvalhais et al., 2014; Schimel et al., 2015) while comprehensively considering various factors from both sides that affect CO₂ in one coordinate frame, which could be a very helpful tool to support policymakers for balancing short-term carbon cycles at regional scales.