Anthropogenic and natural controls on atmospheric δ¹³C-CO₂ variations in the Yangtze River delta: insights from a carbon isotope modeling framework

The atmospheric carbon dioxide (CO₂) mixing ratio and its carbon isotope (δ¹³C-CO₂) composition contain important CO₂ sink and source information spanning from ecosystem to global scales. The observation and simulation for both CO₂ and δ¹³C-CO₂ can be used to constrain regional emissions and better understand the anthropogenic and natural mechanisms that control δ¹³C-CO₂ variations. Such work remains rare for urban environments, especially megacities. Here, we used near-continuous CO₂ and δ¹³C-CO₂ measurements, from September 2013 to August 2015, and inverse modeling to constrain the CO₂ budget and investigate the main factors that dominated δ¹³C-CO₂ variations for the Yangtze River delta (YRD) region, one of the largest anthropogenic CO₂ hotspots and densely populated regions in China. We used the WRF-STILT model framework with category-specified EDGAR v4.3.2 CO₂ inventories to simulate hourly CO₂ mixing ratios and δ¹³C-CO₂, evaluated these simulations with observations, and constrained the total anthropogenic CO₂ emission. We show that (1) top-down and bottom-up estimates of anthropogenic CO₂ emissions agreed well (bias < 6 %) on an annual basis, (2) the WRF-STILT model can generally reproduce the observed diel and seasonal atmospheric δ¹³C-CO₂ variations, and (3) anthropogenic CO₂ emissions played a much larger role than ecosystems in controlling the δ¹³C-CO₂ seasonality. When excluding ecosystem respiration and photosynthetic discrimination in the YRD area, δ¹³C-CO₂ seasonality increased from 1.53 ‰ to 1.66 ‰. (4) Atmospheric transport processes in summer amplified the cement CO₂ enhancement proportions in the YRD area, which dominated monthly δ_s (the mixture of δ¹³C-CO₂ from all regional end-members) variations. These findings show that the combination of long-term atmospheric carbon isotope observations and inverse modeling can provide a powerful constraint on the carbon cycle of these complex megacities.

1 Introduction

Urban landscapes account for 70 % of global CO₂ emissions and represent less than 3 % of Earth's land area (Seto et al., 2014). Such CO₂ hotspots play a dominant role in controlling the rise in atmospheric CO₂ concentrations, which exceeded 412 ppm in December 2019 for global monthly average observations (https://www.esrl.noaa.gov/gmd/ccgg/trends/, last access: 1 August 2020). Furthermore, the carbon isotope ratio of CO₂ (i.e., δ¹³C = ¹³C $/$ ¹²C ratio in delta notation) at the representative Mauna Loa site, USA, has steadily decreased to around −8.5 ‰ in December 2019 (https://www.esrl.noaa.gov/, last access: 1 August 2020). Anthropogenic CO₂ emission is produced from fossil fuel burning and cement production. As the urban population is expected to increase by 2.5 to 6 billion people in 2050, anthropogenic CO₂ emissions are projected to increase dramatically, especially in developing regions and countries (Sargent et al., 2018; Ribeiro et al., 2019). Under such a scenario, the observations of atmospheric CO₂ and δ¹³C-CO₂ in urban landscapes are of great importance to monitoring these potential CO₂ emissions hotspots (Lauvaux et al., 2016; Nathan et al., 2018; Graven et al., 2018; Pillai et al., 2016; Staufer et al., 2016).

Countries are required to report their CO₂ emissions according to the Intergovernmental Panel on Climate Change guidelines (IPCC, 2019), and many bottom-up methods have long been used to estimate CO₂ emissions worldwide, but such methods have high uncertainties for CO₂ emissions at regional (20 %) to city (50 % to 250 %) scales (Gately and Hutyra, 2017; Gately et al., 2015). These large uncertainties are propagated into the estimation of biological fluxes in atmospheric inversions (Zhang et al., 2014; Jiang et al., 2014; Thompson et al., 2016). By using CO₂ observations, the top-down atmospheric inversion approach is a useful tool to evaluate bottom-up inventories (Graven et al., 2018; L. Hu et al., 2019; Lauvaux et al., 2016; Nathan et al., 2018). Previous research has shown that additional information, such as data on atmospheric Δ¹⁴CO₂-CO₂, δ¹³C-CO₂, and CO, is needed to better distinguish CO₂ emissions from different sources and to assess their uncertainties (Chen et al., 2017; Graven et al., 2018; Nathan et al., 2018; Cui et al., 2019). The use of hourly δ¹³C-CO₂ observation in urban areas remains rare in inversion studies, yet such observations contain invaluable information of anthropogenic CO₂ from different categories.

Traditional estimates of δ¹³C-CO₂ using isotope ratio mass spectrometry (IRMS) are very limited because flask air sample collection requires long preparation time and is expensive. Consequently, there is a lack of high temporal and long-term observations of δ¹³C-CO₂ (Sturm et al., 2006). Isotope ratio infrared spectroscopy (IRIS) technology has overcome these limitations. As a result, in situ air sample analyses using IRIS analyzers result in dense time series of δ¹³C-CO₂. However, most of the established long-term IRMS and IRIS δ¹³C-CO₂ measurement sites are representative of “background”, natural, or agricultural ecosystems at locations far away from urban landscapes (Chen et al., 2017; Griffis, 2013).

To date, long-term (> 1-year) and continuous observations of both CO₂ and δ¹³C-CO₂ have been reported for only five cities, including Bern, Switzerland (Sturm et al., 2006), Boston, USA (McManus et al., 2010), Salt Lake City, USA (Pataki et al., 2006), Beijing, China (Pang et al., 2016), and Nanjing, China (Xu et al., 2017). In these previous investigations, significant diel and seasonal variations of δ¹³C-CO₂ have been observed; these patterns were modulated by fossil fuel combustion, plant respiration and photosynthesis, and changes in the height of the atmospheric boundary layer (Sturm et al., 2006; Guha and Ghosh, 2010). No study has quantified the impact of each factor on the seasonal variation of δ¹³C-CO₂. This represents an important knowledge gap in understanding the underlying mechanisms of carbon cycling in complex urban ecosystems.

The traditional δ¹³C-CO₂ isotope partitioning methods (including the Miller–Tans and Keeling plot approaches) have been used to constrain different CO₂ sources worldwide (Keeling, 1960; Vardag et al., 2015; Newman et al., 2016; Pang et al., 2016; Xu et al., 2017). These methods are based on the assumption that partitioned atmospheric CO₂ enhancement components from different sources can represent CO₂ emissions in the “target area” (Miller et al., 2003; Ballantyne et al., 2011). Carbon dioxide emissions are highly inhomogeneous at the urban scale, with extremely strong point/line sources, and the final partitioning results are highly uncertain without considerations of source footprint characteristics (Gately and Hutyra, 2017; Cui et al., 2019; Martin et al., 2019). Atmospheric transport models can help to resolve such problems, and the coupling of atmospheric transport models with isotope observations has recently been applied in global and regional CO₂ partitioning studies (Chen et al., 2017; Cui et al., 2019; Graven et al., 2018; Hu et al., 2018b). Although urban CO₂ inversions have been applied successfully in several studies in Europe and the United States (Bréon et al., 2015; Turnbull et al., 2015; Pillai et al., 2016; Brioude et al., 2013; Turner et al., 2016), urban CO₂ inversions in China are rare (Berezin et al., 2013; Hu et al., 2018a; Worden et al., 2012), presumably because of the scarcity of high-quality δ¹³C-CO₂ and CO₂ observations.

The Yangtze River delta (YRD) ranks as one of the most densely populated regions in the world and is an important anthropogenic CO₂ hotspot. Major anthropogenic sources include the power industry, oil refineries/transformation and cement production. Having the largest source of cement-derived CO₂ production across China and the world (Cai et al., 2015), the YRD contributed 20 % of national cement production, nearly 12 % of the world's total cement output in 2014 (USGS, 2014; Xu et al., 2017; Yang et al., 2017). In addition to anthropogenic factors, natural ecosystems and croplands act as significant CO₂ sinks and sources within the YRD. Independent quantification of the fossil and cement CO₂ emission and assessment of their impact on atmospheric δ¹³C-CO₂ have the potential to improve our understanding of urban CO₂ cycling. Further, the observations and simulations of both atmospheric CO₂ and δ¹³C-CO₂ can help us relate atmospheric CO₂ dynamics to future emission control strategies.

Here, we combine long-term (>2-year) CO₂ and δ¹³C-CO₂ observations with atmospheric transport model simulations to study urban atmospheric CO₂ and δ¹³C-CO₂ variations. The objectives were to (1) constrain anthropogenic CO₂ emissions and determine the main sources of uncertainty for δ¹³C-CO₂ simulations and (2) quantify the relative contributions of each factor (i.e., background, anthropogenic CO₂ emissions especially for cement production, ecosystem photosynthesis and respiration) to seasonal variations of atmospheric δ¹³C-CO₂.

2 Materials and methods

2.1 Observations of atmospheric CO₂ mixing ratio, δ¹³C-CO₂ and supporting variables

The observation site is located on the Nanjing University of Information Science and Technology campus (hereafter NUIST, 32^∘12^′ N, 118^∘43^′ E; green dot in Fig. 1a). Continuous atmospheric CO₂ mixing ratios and δ¹³C-CO₂ were measured at a height of 34 m above ground with an IRIS analyzer (model G1101-i, Picarro Inc., Sunnyvale, CA). The observation period extended from September 2013 to August 2015. Calibrations for CO₂ mixing ratio and δ¹³C-CO₂ were conducted with standard gases traceable to NOAA/GML (NOAA Global Monitoring Laboratory) standards. Calibration details are provided by Xu et al. (2017). Based on Allan variance analyses, the hourly precisions of CO₂ and δ¹³C-CO₂ were 0.07 ppm and 0.05 ‰, respectively. We note that the δ¹³C-CO₂ IRIS (model G1101-i) measurements are sensitive to water vapor concentration. Sensitivity tests reveal that the δ¹³C-CO₂ IRIS measurements are biased high (less than 0.74 ‰) when water vapor mole fraction exceeds 2 %. The data presented here have been corrected following the procedures outlined in Xu et al. (2017).

Figure 1 (a) Weather Research and Forecasting Model simulation domains and the location of the WLG site; the different region colors represent three domains. (b) Cement production distribution in the YRD and eastern China. Both the green dot in (a) and the red star in (b) are the NUIST observation site.

We separated the 2-year study period into seasons (autumn: September, October, November; winter: December, January, February; spring: March, April, May; summer: June, July, August). Further, for an annual comparison, we examined the period from September 2013 to August 2014 (year 2014) versus September 2014 to August 2015 (year 2015).

The YRD is a cement production hotspot in China (Fig. 1b). It had a total population of 190 million in 2018 (Fig. 2a), with 24.2 million in the city of Shanghai, 9.8 million in Hangzhou (provincial capital of Zhejiang), 8.4 million in Nanjing (provincial capital of Jiangsu), and 8.1 million in Hefei (provincial capital of Anhui). The CO₂-related production data (i.e., cement) and energy consumption data (i.e., coal and natural gas) were obtained from local official sources using the same method described in Shen et al. (2014).

To examine the effects of plant photosynthesis on atmospheric CO₂ variations, we used the NDVI (Normalized Difference Vegetation Index), SIF (solar-induced chlorophyll fluorescence) and GPP (gross primary productivity) information. These three products have a global distribution with a spatial resolution of 0.05^∘ by 0.05^∘. The NDVI has a temporal resolution of 16 d, and SIF and GPP products have a temporal resolution of 8 d (Li and Xiao, 2019; http://globalecology.unh.edu/data/, last access: 5 March 2020). Land-use and land-cover classification in the Yangtze River delta for 2014 was applied by using NDVI data from MOD13A2.

2.2 Simulation of atmospheric δ¹³C-CO₂

2.2.1 General equations

The simulation of atmospheric δ¹³C-CO₂ is based on mass conservation. First, we briefly describe the simulation of atmospheric CO₂ mixing ratios (more details are provided in Sect. 2.2.2), following the previous work of Hu et al. (2018b), where atmospheric CO₂ was simulated (CO_{2_sim}) as the sum of background (CO_{2_bg}) and the contribution from all regional sources/sinks ([ΔCO_{2_sim}]_i), as

$$\begin{array}{}\text{(1)}& {\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}={\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{bg}}+\sum _{i=\mathrm{1}}^n[\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}{]}_i.\end{array}$$

Note that ΔCO₂ is the sum of all simulated sources/sinks [ΔCO_{2_sim}]_i and represents the total simulated CO₂ enhancement. We use ΔCO_{2_obs} as the observed CO₂ total enhancement, which can be calculated by using the CO₂ observation minus the CO₂ background values. Based on mass conservation, we estimated the ¹³CO₂ composition by multiplying the left- and right-hand sides of Eq. (1) by δ¹³C:

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}\mathrm{\_}\mathrm{sim}}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{bg}}\times {\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{bg}}+\sum _{i=\mathrm{1}}^n{\mathit{\delta }}_i^{\mathrm{13}}\times [\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}{]}_i}{{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}}},\end{array}$$

where δ¹³C_{a_sim} and δ¹³C_bg represent the simulated atmospheric δ¹³C-CO₂ and background δ¹³CO₂, and ${\mathit{\delta }}_i^{\mathrm{13}}$ is the δ¹³C-CO₂ for end-member i (including anthropogenic and biological source categories). The δ¹³C-CO₂ contributions from all regional sources/sinks can be further reformatted as Eq. (3):

$$\begin{array}{}\text{(3)}& \sum _{i=\mathrm{1}}^n{\mathit{\delta }}_i^{\mathrm{13}}\times [\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}{]}_i={\mathit{\delta }}_{\mathrm{s}\mathrm{\_}\mathrm{sim}}\times \sum _{i=\mathrm{1}}^n[\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}{]}_i,\end{array}$$

where δ_{s_sim} is the simulated enhancement-weighted mean of all regional end-members. We use δ_s as the observed term to distinguish it from δ_{s_sim} (Newman et al., 2008), which will be described in detail in Sect. 2.2.5. The product on the right-hand side of Eq. (3) is the simulated regional source term that is added to the background value and contains both enhancement and δ¹³C-CO₂ signals contributed by different CO₂ sources/sinks. This product can also be treated as an observed term when using the derived δ_{s_obs} and observed δCO_{2_obs} values.

To date, there are no available global δ¹³C-CO₂ background products, and the choice of δ¹³C_bg is essential for simulating δ¹³C_a. Here, we apply three strategies. First, we used discrete δ¹³C-CO₂ flask observations at Mount Waliguan (hereafter WLG, 36^∘17^′ N, 100^∘54^′ E; https://www.esrl.noaa.gov/gmd/dv/data/, last access: 31 December 2019) to represent the δ¹³C-CO₂ background signal at our site. These observations were measured at weekly intervals to the end of 2015. A digital filtering curve-fitting (CCGCRV) regression method was applied to derive hourly background values following Thoning et al. (1989). There are, however, reasons why WLG may not be an ideal background site for our study domain. For example, based on the previous simulation results for the CO₂ background sources, most of the back trajectories originate from the free atmosphere or 1000 m higher above the ground (C. Hu et al., 2019). Further, the footprint at the northern/western edge of domain 1 is relatively small, indicating that most back trajectories were observed above the planetary boundary layer height (hereafter PBLH). Here, the WLG observations were made near the surface. Further, WLG is not located at the border of our simulation domain 1. Therefore, the strong vertical δ¹³C-CO₂ gradients between the boundary layer and the free tropospheric atmosphere (Chen et al., 2006; Guha and Gosh, 2010; Sturm et al., 2013) can cause a low bias in the δ¹³C-CO₂ background when using this approach.

In the second approach, the δ¹³C-CO₂ background signal was estimated with wintertime “clean” air CO₂ and δ¹³C-CO₂ observations at the NUIST site, using the following equation:

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{bg}}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{a}}\times {\mathrm{CO}}_{\mathrm{2}}-\sum _{i=\mathrm{1}}^n{\mathit{\delta }}_i^{\mathrm{13}}\times [\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}{]}_i}{{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{bg}}}},\end{array}$$

where δ¹³C_a and CO₂ represent atmospheric δ¹³C-CO₂ and CO₂ observations at the NUIST site under clean conditions. Note that δ¹³C_a represents the observed δ¹³C-CO₂, not the simulated δ¹³C-CO₂ (δ¹³C_{a_sim}) as shown in Eq. (2). [ΔCO_{2_sim}]_i is the simulated category-specified CO₂ enhancements. We defined clean conditions as the bottom 5 % wintertime CO₂ observations to minimize simulated CO₂ enhancement errors from both biological and anthropogenic CO₂ simulations on δ¹³C-CO₂ background calculation. The CO_{2_bg} is obtained from heights 1000 m above ground level (see Sect. 2.2.3).

In the third approach, we avoid the use of modeled [ΔCO_{2_sim}]_i results and replaced the simulated regional source term in Eq. 4 with observed δ_{s_obs}×CO_{2_obs}, as described in Eq. (3), and used the Miller–Tans regression method to calculate monthly δ_{s_obs}. This approach does not require simulation of [ΔCO₂]_i or the corresponding δ¹³C-CO₂ signals. The hourly δ¹³C-CO₂ background value can be derived by using δ_{s_obs}, CO₂ background, observed atmospheric δ¹³C_a and CO₂ (see details in Sect. 2.3 and the Supplement). Comparison of these three strategies will be evaluated and discussed in Sect. 3.2.1. Similar methods used to derive other background tracers have included CO₂ (Alden et al., 2016; Verhulst et al., 2017), CO (Wang et al., 2010; Ruckstuhl et al., 2012) and CH₄ (Zhao et al., 2009; Verhulst et al., 2017; C. Hu et al., 2019). To analyze the controlling factors for the δ¹³C-CO₂ seasonality, the CCGCRV (a digital filtering curve-fitting program developed by the Carbon Cycle Group, NOAA, USA) regression was applied to the background, observations, and simulations. Finally, we derived CCGCRV curve-fitting lines by using 11 regressed parameters, which were based on the hourly time series of observations/simulations, and defined the difference between peak and trough in 1 year as the seasonality of δ¹³C-CO₂.

2.2.2 Simulation of atmospheric CO₂ mixing ratios

In Eq. (1), the CO_{2_bg} is obtained from the Carbon Tracker 2016 product, which provides global CO₂ distributions from the ground level up to a height of 50 km. We used the averaged concentration above the latitude and longitude where the released particles entered study domain 1 (Fig. 1a). The variable ΔCO_{2_sim} was derived by multiplying the simulated hourly footprint function by the hourly CO₂ fluxes (Hu et al., 2018a, b). Considering the diurnal variations of both anthropogenic and biological CO₂ fluxes, 168 footprints were obtained representing each simulated hour. This accounted for the back trajectory of particle movement for 168 h (i.e., 24 h per day for 7 d) of transport. The 168 footprints are multiplied by the corresponding hourly CO₂ flux. The CO₂ fluxes contain anthropogenic CO₂ emissions, biological CO₂ flux and biomass burning. Here the anthropogenic CO₂ emission sources include the power industry, combustion for manufacturing, non-metallic mineral production (cement), oil refineries/transformation industry, energy for building and road transportation. Theoretically, ΔCO_{2_sim} represents the CO₂ changes contributed by every pixel within the simulated domain. As shown by Hu et al. (2018a), most of the ΔCO_{2_sim} is contributed by sink/source activity within the YRD area. In order to quantify the relative contributions within the YRD area, we separated the study domain into five zones based on provincial administrative boundaries including Jiangsu, Anhui, Zhejiang, Shanghai, and the remaining area outside the YRD (Fig. 2). The modeled CO₂ was calculated as follows:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{sim}}=\sum _{i=\mathrm{1}}^{\mathrm{168}}{\mathrm{flux}}_i\times {\mathrm{footprint}}_i,\end{array}$$

where flux_i (units: mol m⁻² s⁻¹) corresponds to each CO₂ flux category simulated for each domain for a specific hour i, and footprint (units: ppm m² s µmol⁻¹) is the model-simulated sensitivity of observed CO₂ enhancement to flux changes in each pixel. The i contains the hourly footprint during the trajectory of particle movement for 168 h as described above. The CO₂ enhancements from each of the five zones were simulated by multiplying CO₂ emissions in each province by the corresponding footprint.

Figure 2 (a) Annual anthropogenic CO₂ emissions for the study domain (units: mol m⁻² s⁻¹) and population density in four megacities (units: people per hectare) including Nanjing, Hefei, Zhejiang, and Shanghai for the year of 2015 and (b) the 2-year average concentration footprint (units: ppm m² s mol⁻¹).

2.2.3 WRF-STILT model configuration

The Stochastic Time-Inverted Lagrangian Transport (hereafter STILT) model was used to generate the above footprint, which is defined as the sensitivity of atmospheric CO₂ enhancement to the upwind flux at the receptor site (observation site). The meteorological fields used to drive the STILT model were simulated with the Weather Research and Forecasting Model (WRF3.5) at high spatial and temporal resolutions. The innermost nested domain (D3, 3 km × 3 km, Fig. 1) contains the YRD area, where the most sensitive footprint is located, and the intermediate domain (D2, 9 km × 9 km) and outermost domain (D1, 27 km × 27 km) represent eastern China and central and eastern China, respectively. The same physical schemes and parameter setup for the WRF meteorological field simulation and the domain in the STILT model have been used previously for inverse analyses (C. Hu et al., 2019). These previous studies at the NUIST observation site have shown very good performance in simulating the meteorological fields, which is essential for reliable STILT simulations. The hourly footprint was simulated by releasing 500 particles from the NUIST measurement site and tracking their backward locations every 5 min for a period of 7 d. Particle numbers and their residence time within half of the PBLH were used to calculate the footprint over the 7 d period. For the CO₂ background of each hour, we tracked the sources of air particles' back trajectory for 7 d and defined these CO₂ mixing ratios in Carbon Tracker as the hourly CO₂ background values (Peters et al., 2007).

2.2.4 A priori anthropogenic CO₂ emissions and net ecosystem exchange

The Emission Database for Global Atmospheric Research (EDGAR v4.3.2) inventory was selected as the a priori anthropogenic CO₂ emissions (Fig. 2a), which is based on the International Energy Agency's (IEA's) energy budget statistics and provides detailed CO₂ source maps (29 categories, including both organic and fossil emissions, IEA, 2012) with global coverage at high spatial resolution (0.1^∘ × 0.1^∘). The EDGAR CO₂ emissions are the most up-to-date global inventory with sectoral detail (Janssens-Maenhout et al., 2017; Schneising et al., 2013). Other inventories, including the Fossil Fuel Data Assimilation System (FFDAS, Rayner et al., 2010) and the Open-source Data Inventory for Anthropogenic CO₂ (ODIAC, Oda et al., 2018), also provide global CO₂ emissions. However, these inventories only provide total CO₂ emissions or have very limited emission categories, which limits our ability to provide isotope end-member information. EDGAR v4.3.2 provides emission estimates at a monthly timescale. Here, we applied hourly scaling factors for different categories following Hu et al. (2018a). EDGAR v4.3.2 with monthly resolution is available only for 2010. We assume that each CO₂ category changes linearly from its 2010 value (Peters et al., 2007) and apply an annual scaling factor of 1.145 to derive CO₂ emissions for 2014 and 2015. This scaling factor is based on Carbon Tracker, dividing the same anthropogenic CO₂ emissions for the YRD in years 2014–2015 by that in 2010.

The biological flux or net ecosystem CO₂ exchange (NEE) and biomass burning CO₂ emissions come from Carbon Tracker a posteriori flux at 3 h intervals and at a spatial resolution of 1^∘ × 1^∘. Because NEE is much smaller than the anthropogenic CO₂ emissions in such densely developed urban landscapes, we homogeneously distributed this flux at a spatial resolution of 0.1^∘ within each grid to match the footprint.

2.2.5 Simulation of the carbon isotope ratio of all sources (δ_{s_sim})

The carbon isotope ratio of all the surface sources was calculated as (Newman et al., 2008)

$$\begin{array}{}\text{(6)}& \sum _{i=\mathrm{1}}^n{\mathit{\delta }}_i\times p_i={\mathit{\delta }}_{\mathrm{s}\mathrm{\_}\mathrm{sim}},\end{array}$$

where δ_i is the δ¹³C-CO₂ value from source category i, and p_i is the corresponding enhancement proportion (i.e., proportions of a specific enhancement i to total CO₂ enhancement). We define δ_{s_sim} as the simulated carbon isotope ratio of all sources to differentiate it from the observed δ_{s_obs}. Based on fossil fuel usage characteristics in the YRD, we reassigned the EDGAR v4.3.2 categories according to fuel types. Coal was the fuel type for manufacturing, oil for oil refinery, natural gas for buildings, and diesel and gasoline for transportation. The power industry consumed 5 % natural gas and 95 % coal based on local activity data in the YRD (State Statistical Bureau, 2016). The non-metallic mineral production was mainly for cement. Since there is a lack of detailed information for non-metallic mineral production, we simply attributed 100 % of it to cement production. Chemical processes were mainly ammonia synthesis. Based on a literature review and our previous work (Xu et al., 2017), typical δ¹³C-CO₂ values for natural gas (−39.06 ‰ ± 1.07 ‰), coal (−25.46 ‰ ± 0.39 ‰), fuel oil (−29.32 ‰ ± 0.15 ‰), gasoline (−28.69 ‰ ± 0.50 ‰), ammonia synthesis (−28.18 ‰ ± 0.55 ‰), diesel (−28.93 ‰ ± 0.26 ‰), pig iron (−24.90 ‰ ± 0.40 ‰), crude steel (−25.28 ‰ ± 0.40 ‰), cement (0 ‰ ± 0.30 ‰), and biofuel combustion and biological emissions (−28.20 ‰ ± 1.00 ‰) were used in this study. We also applied a value of −28.20 ‰ for photosynthesis (Griffis et al., 2008; Lai et al., 2014) because the YRD is a region dominated by C₃ plants. Since CO₂ emissions associated with human respiration (Prairie and Duarte, 2007; Turnbull et al., 2015; Miller et al., 2020) are relatively small (3.7 % of anthropogenic emissions in the YRD area, Xu et al., 2017) and given that the local food diet is dominated by C₃ grains that have a similar δ¹³C-CO₂ value to the biological CO₂ flux of −28.20 ‰, we assume it has the same isotope signals as local C₃ plants and ecosystem respiration. Further, the biological CO₂ flux from the Carbon Tracker assimilation system considered anthropogenic emissions to be fixed and attributed the remainder to the biological CO₂ flux (Peters et al., 2007). Consequently, we believe the uncertainty in the biological CO₂ flux will include the small proportion of human respiration.

To evaluate the simulated δ_{s_sim}, we applied the Miller–Tans and Keeling plot approaches to derive δ_{s_obs} from the observed concentration and atmospheric ¹³CO₂-CO₂ (Xu et al., 2017). We then used the results to evaluate the calculations made with Eq. (6).

2.3 Independent IPCC method for anthropogenic CO₂ emissions

Large differences among inventories have been previously found even for the same region (Berezin et al., 2013; Andrew, 2018). For comparison with the EDGAR v4.3.2 inventory results, we derived the anthropogenic CO₂ emissions by using an independent IPCC method. Here, we illustrate the calculation for cement CO₂ emissions. Note that the IPCC only recommended an EF for clinker, which is an intermediate product of cement. To calculate cement CO₂ emissions, we need to calculate it based on clinker production, as shown in Eq. (7):

$$\begin{array}{}\text{(7)}& {\mathrm{CO}}_{\mathrm{2}}\left[\mathrm{cement}\right]=M_{\mathrm{cement}}\times C_{\mathrm{clinker}}\times {\mathrm{EF}}_{\mathrm{clinker}},\end{array}$$

where CO₂ [cement] is the chemical process CO₂ emissions for cement production, M_cement is the production of cement, C_clinker represents the clinker-to-cement ratio (%), and EF_clinker is the CO₂ emission factor for clinker production. The IPCC recommended an EF_clinker value of 0.52±0.01 t CO₂ per tonne clinker produced, where CaO content for clinker is assumed to be 65 % with 100 % CaO from calcium carbonate material (IPCC, 2019). The EF appears to be well constrained, showing little variation among provinces with mean values ranging from 0.512 to 0.525 (Yang et al., 2017). For the C_clinker values, it generally showed a decreasing trend from 64.5 % in 2004 to 56.9 % in 2015 for all of China (Fig. S1 in the Supplement), with an average value of 57.0 % during 2014 and 2015.

2.4 Multiplicative scaling factor method

To quantify anthropogenic CO₂ emissions and to compare them with EDGAR products, we first derived the monthly scaling factors for anthropogenic CO₂ emissions using a multiplicative scaling factor (hereafter MSF) method (Sargent et al., 2018; He et al., 2020) and then obtained annual averages. The monthly scaling factors (SFs) were calculated as

$$\begin{array}{}\text{(8)}& \mathrm{MSF}={\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{obs}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{bg}}-\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{bio}}-\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{fire}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}\mathrm{\_}\mathrm{anthro}}}},\end{array}$$

where CO_{2_obs}, ΔCO_{2_bio}, ΔCO_{2_fire} and ΔCO_{2_anthro} represent observed CO₂ mixing ratios, simulated CO₂ enhancements contributed by biological flux, biomass burning, and anthropogenic emissions, respectively. Uncertainties of all factors on the final MSFs were calculated based on Monte Carlo methods, where the normal sample probability distribution was applied and the upper 97.5 % and lower 2.5 % of the values were considered to be the uncertainty for MSF (Cao et al., 2016).

3 Results and discussion

3.1 Evaluation of hourly CO₂ mixing ratios

3.1.1 Hourly and monthly CO₂ mixing ratio comparisons

This section examines the general performance of simulating hourly CO₂ mixing ratios. The 2-year average hourly footprint is shown in Fig. 2b, where the source area (blue–red) indicates strong sensitivity of the CO₂ observations to regional sources. This footprint shape is representative of the YRD area. To quantify the relative contributions from each province, we calculated CO₂ enhancements contributed by Anhui, Jiangsu, Zhejiang, Shanghai, and the remaining area outside of the YRD, respectively. The results indicate that Jiangsu contributed approximately 80 % of the total enhancement (discussed further in Sect. 3.1.2). Comparisons between simulated and observed hourly CO₂ mixing ratios are displayed in Fig. 3a for both years. For all hourly data in each year, the model versus observation correlation coefficient (R) was R=0.38 (n=8204, P<0.001) and RMSE = 29.44 ppm for 2014 and R= 0.35 (n=7262, P<0.001) and RMSE = 30.22 ppm for 2015. These results indicate that the model can simulate the synoptic and diel CO₂ variations over the 2-year period. The model also captured the monthly and seasonal variations of CO₂ mixing ratios (daily averages are shown in Fig. S2). The simulations captured the trend of rising CO₂ mixing ratios after October and the drawdown of CO₂ to the background value during the summer.

Figure 3 (a) Comparisons of hourly CO₂ mixing ratios between observations and model simulation from September 2013 to August 2015 and monthly averages for (b) whole day, (c) nighttime (22:00–06:00, local time) and (d) daytime (10:00–16:00). Model results (red), observations (black), and background (grey).

Figure 3b–d illustrate the average monthly daily, nighttime (22:00–06:00, local time), and daytime (10:00–16:00) CO₂ mixing ratios. These monthly values contain the effects of atmospheric transport, background and variations in CO₂ emissions. The observed and simulated CO₂ mixing ratios showed a significant increase from September 2013 to January 2014. Here, the CO₂ mixing ratios increased by 16.0 ppm according to the model results and 17.2 ppm according to the observations. The background values increased by 8.1 ppm and accounted for 47 % of the total CO₂ increase, and the net CO₂ flux (a priori) for the YRD increased by 15 %. We attributed the remaining 38 % increase to changes in atmospheric transport processes including lower PBLH in January 2014 than in September 2013. To quantify how variations in PBLH affected CO₂ mixing ratios, we compared the simulated monthly anthropogenic CO₂ enhancement differences in the same months of different years to eliminate the influence of monthly emission variations on CO₂ enhancements. Twelve monthly paired values were used and are shown in Fig. 4. This analysis indicates that atmospheric CO₂ mixing ratios decreased by about 3.7 ppm for an increase in PBLH by 100 m. We also note that there were 2 months (March and August) that fall far below this trend, implying that changes in the monthly footprints (source area) can also play an important role.

Figure 4 Relation between monthly PBLH and change in the CO₂ mixing ratio; here, these dots represent the difference of monthly averages in 2 different years for all hours.

On an annual timescale, the simulated average CO₂ mixing ratios were 436.63 and 437.11 ppm for 2014 and 2015, respectively. Since the anthropogenic CO₂ emissions used in the model are the same for both years, the simulated annual average CO₂ difference can be used to quantify the influence associated with meteorological factors and ecosystem carbon cycling. Between these 2 years, the CO₂ background increased by 1.78 ppm, and the biological enhancement decreased by 1.04 ppm from 2014 to 2015. The remaining 0.26 ppm change between 2014 and 2015 indicates a relatively small meteorological effect for the annual averages, such as a slight change in the dominant wind direction or a PBLH difference.

The simulated annual average NEE CO₂ enhancements were 2.64 and 1.60 ppm for the respective years. For comparison, the annual average anthropogenic enhancements were 36.20 and 34.90 ppm for 2014 and 2015, respectively. The monthly NEE enhancement varied from −0.1 ppm in May 2015 to +6.0 ppm in July 2014, indicating NEE contributes positively for enhancement in most months (Fig. 5a), even though the sign of the monthly averaged NEE flux in summer was negative (sinks). This positive contribution was mainly caused by diel PBLH variations between daytime (smaller negative enhancement) and nighttime (larger positive enhancement). To further evaluate the impact of plant photosynthetic activity on the regional CO₂ cycle, we examined the NDVI, SIF and GPP seasonal patterns (Fig. 5d–f). These three datasets revealed two peaks during each year, which is related to increased photosynthetic activity. The first peak occurred in May and the second in August–September, corresponding to the growing season of wheat and corn/rice, respectively (Deng et al., 2015). We note that GPP was derived from SIF, and as a result, they share a similar seasonal cycle. The land-use classification in the YRD for 2014 (Fig. S3) shows that the northern YRD is dominated by agricultural land and the south dominated by forest land, and our observation site was more surrounded by agricultural land, which corresponded well to observed NDVI, SIF and GPP seasonal patterns. The peak SIF and GPP signals during the summer were about 20 times greater than during the winter. Consequently, we can ignore the potential influence of photosynthetic activity on the regional CO₂ enhancements during the non-growing seasons.

Figure 5 (a) Comparisons of simulated and observed CO₂ enhancement; note that “model” represents the sum of both anthropogenic and biological CO₂ enhancement simulations, (b) CO₂ enhancement contributions from different provinces, and (c) simulated anthropogenic CO₂ enhancement proportion for the main sources. Time series (2013 to 2015) of (d) NDVI, (e) SIF, and (f) GPP. The distance indicates the radius of the area centered with the NUIST observation site, and the NDVI, SIF, and GPP values are averages in these areas.

3.1.2 Components of urban CO₂ enhancement

Here, we diagnose the source contributions to the urban CO₂ enhancement. The observed anthropogenic CO₂ enhancements, which were derived by subtracting CO₂ background and simulated biological enhancement from CO₂ concentration observations, were 38.36 ± 3.32 and 37.89 ± 2.80 ppm for 2014 and 2015, respectively. Here, the uncertainty of the observed anthropogenic CO₂ enhancements was calculated by prescribing a 2 ppm potential bias for the Carbon Tracker CO₂ fields and 50 % to the simulated biological CO₂ enhancement (Hu et al., 2018b). The corresponding simulated anthropogenic CO₂ enhancements were 36.20 and 34.90 ppm. In comparison with the simulated biological CO₂ enhancements displayed in Fig. 5a, both the observed and simulated CO₂ enhancements are indicative of a large anthropogenic (fossil fuel and cement production) CO₂ emission from the YRD.

Previous studies have also investigated urban CO₂ enhancements from a relatively broad range of developed environments worldwide. Verhulst et al. (2017) measured CO₂ mixing ratios at seven sites in Los Angeles, USA, and concluded that the mean annual enhancement varied between 2.0 and 30.8 ppm, which is considerably lower than our findings. Another study in Washington DC, USA, in February and July 2013 showed that the CO₂ enhancement was less than 20 ppm (Mueller et al., 2018). The urban CO₂ observations and modeling study by Martin et al. (2019) at three urban sites in the eastern USA showed an enhancement of ∼ 21 ppm in February 2013, substantially lower (by ∼ 20 ppm) than our observations. The measurements at an urban–industrial complex site in Rotterdam, Netherlands, indicated a CO₂ enhancement of only 11 ppm for October to December 2014 (Super et al., 2017). Our enhancements were significantly higher than all of these previous reports of other urban areas.

The anthropogenic components and source area contributions are displayed in Fig. 5b–c. During the study period the average anthropogenic enhancements were 5.1 %, 80.2 %, 1.9 %, 4.4 %, and 8.5 % for Anhui, Jiangsu, Zhejiang, Shanghai, and the remaining area outside the YRD, respectively. Although Shanghai's area is the smallest within the YRD region and relatively distant (∼ 300 km) from our observation site, its maximum source contribution at times exceeded 50 % (i.e., on 19 September 2013, not shown) via long-distance transport. In general, the power industry, manufacturing, non-metallic mineral production, oil refinery, and other source categories contributed 41.0 %, 21.9 %, 9.3 %, 11.5 %, and 16.3 % to the total anthropogenic CO₂ enhancement, respectively. The proportions of corresponding CO₂ emission categories to the total anthropogenic emissions of the YRD were 39.8 %, 28.4 %, 7.4 %, 4.1 %, and 24.4 %, respectively. The comparisons between the proportions of simulated enhancement and proportions of corresponding CO₂ emissions can illustrate whether CO₂ enhancement partitions are a good tracer for emissions in a complex urban area. We found a relatively large difference between the enhancement proportion and the emission proportion for oil refineries (from 11.5 % to 4.1 %) as compared to other categories. This may be because the power industry, manufacturing and non-metallic mineral production were more homogeneously distributed compared to oil refineries, which were closer to our CO₂ observation site. Further, changes in source footprint caused by wind direction variations likely played an important role.

3.1.3 Constraints on monthly anthropogenic CO₂ emissions

To provide a robust comparison of bottom-up CO₂ emissions for the YRD, we calculated anthropogenic CO₂ emissions from both EDGAR v4.3.2 and with activity data provided by local governments (Table 1) and the default IPCC emission factors (https://www.ipcc-nggip.iges.or.jp/EFDB/, last access: 13 September 2019). The total anthropogenic CO₂ emissions in 2014–2015 were 24.4 × 10¹¹ and 23.5 × 10¹¹ kg according to our own inventory and EDGAR v4.3.2 CO₂, respectively, indicating excellent agreement (within 4 %) between these approaches. We constrained the monthly anthropogenic CO₂ emissions by using the MSF method (Eq. 8) and computed the 12-month average to represent the years of 2014 and 2015. The a posteriori results indicate that the annual scaling factors were 1.03 ± 0.10 for 2014 and 1.06 ± 0.09 for 2015. The monthly scaling factors derived from using daytime and all-day observations are also shown in Fig. S4. These factors vary seasonally, with higher values observed in summer. When using daytime values only, the scaling factors were much larger than the all-day values. This can be seen in Fig. 3 by comparing the simulated and observed CO₂ mixing ratios. We should note here that the larger scaling factors based on the daytime data could be caused by bias in the a priori daily scaling factors used to generate the hourly CO₂ emissions (Hu et al., 2018b), the monthly anthropogenic averages, and bias in negative biological CO₂ enhancement. Since our study is mainly focused on the seasonality of all-day observations, the monthly scaling factors derived from the all-day approach will be used for the following analyses. The anthropogenic CO₂ emissions in year 2015 did not show a significant change compared to 2014, and the overall estimates were within the uncertainty of the estimates. After applying the average scaling factors for 2014 and 2015, the a posteriori anthropogenic CO₂ emissions were 24.6 (± 2.4) × 10¹¹ kg for the YRD area. The application of the MSF method provides an overall constraint on the anthropogenic CO₂ emissions (also displayed in Table 1).

Table 1 Comparisons of cement and all anthropogenic CO₂ emissions among the different methods. “NA” means not available.

The main uncertainties associated with the simulation of hourly CO₂ and δ¹³C-CO₂ are uncertainty in meteorological fields, transport model (i.e., number of released particles), and a priori CO₂ fluxes. At the annual scale the main uncertainty is attributed to the PBLH simulations and a priori anthropogenic CO₂ emissions. The anthropogenic CO₂ emissions biases were < 6 % as described above, and the bias associated with PBLH uncertainty was typically < 13 % (Hu et al., 2018a, b). There, we attribute a 20 % uncertainty to the simulated CO₂ and δ¹³C-CO₂ signals on an annual timescale.

3.2 Simulation of atmospheric δ¹³C-CO₂

3.2.1 Background atmospheric δ¹³C-CO₂

To obtain the best representative δ¹³C-CO₂ background value for the study domain, we examined the values from the three strategies described above (Fig. 6). We also compared the δ¹³C-CO₂ at the WLG background site with observations at NUIST during winters (Fig. S5). This was performed to help simplify the comparison by removing the effects of plant photosynthetic discrimination. The δ¹³C-CO₂ at the WLG site was relatively more depleted in the heavy carbon isotope (or negative, by up to 0.5 ‰) than that observed at NUIST for many periods. Theoretically, there are two key factors that can cause the urban atmospheric δ¹³C-CO₂ to be relatively more enriched in the heavy carbon isotope (or positive) compared to the background values, including (1) discrimination associated with ecosystem photosynthesis and (2) enrichment of the isotopic signature associated with the CO₂ derived from cement production. As shown earlier, the biological CO₂ enhancement was positive in winter, which implies a positive biological CO₂ signal where ecosystem respiration is more important than photosynthesis. Further, sensitivity tests for cement CO₂ sources showed its influence is much smaller than the observed difference in Fig. S5 (discussed in Sect. 3.3.3). Based on the above analyses and methods introduced in Sect. 2.3, we concluded that the WLG δ¹³C-CO₂ signal is not an ideal choice for representing the background value. The wintertime δ¹³C-CO₂ background values, based on strategy 2, were −7.78 ‰ and −7.61 ‰ for 2013–2014 and 2014–2015, respectively (Fig. 6). The corresponding values, based on strategy 3, were −7.70 ‰ and −7.53 ‰. These background values are more enriched compared to the WLG observations by 0.80 ‰ to 1.01 ‰. These derived values agree well with the monthly δ¹³C-CO₂ simulation results of Chen et al. (2006), who showed that δ¹³C-CO₂ is 0.6 ‰ higher above the PBL than in the surface layer near the ground. Recently, Ghasemifard et al. (2019) showed that hourly δ¹³C-CO₂ values at the Zugspitze, the highest (2650 m) mountain in Germany, varied between −7 ‰ and −12 ‰ in the winter for 2013. During two especially clean air events (in October and February) at the Zugspitze, the δ¹³C-CO₂ was approximately −7 ‰, during which the CO₂ mixing ratios varied between 390 and 395 ppm. This is consistent with our estimates using strategies 2 and 3. Based on the evidence presented above, we believe that strategy 3 is the most robust way to derive a background δ¹³C-CO₂ for the study domain.

Figure 6 Comparisons among three strategies for calculating the background δ¹³C-CO₂. Strategy 1 (WLG discrete: weekly discrete observations at the WLG site, WLG CCGCRV: derived hourly data with WLG observations and CCGCRV method), strategy 2 (calculated by choosing clean air in winter), and strategy 3 (M-T method: derived results with observations and M-T approach, M-T CCGCRV: derived hourly results with the M-T approach and CCGCRV method; see details in Sect. 2.2.1).

3.2.2 Evaluation of δ¹³C-CO₂ simulations

Figure 7a shows the hourly δ¹³C-CO₂ simulations over a 2-year period. To the best of our knowledge, this is the first time that δ¹³C-CO₂ has been simulated at an hourly timescale for an urban region. The simulations are consistent with the observations at daily, monthly and annual timescales, where the average values of observations (simulations) were −8.69 ‰ (−8.68 ‰) and −8.52 ‰ (−8.45 ‰) for 2014 and 2015, respectively. The corresponding correlations were R=0.54 (P<0.001) and R=0.52 (P<0.001). The root mean square error between observations and simulations was 1.07 ‰ for 2014 and 1.10 ‰ for 2015 (Table 2). Further, the observed and simulated δ¹³C-CO₂ values showed seasonal variations that increased in summer and decreased in winter. This pattern mirrored the CO₂ mixing ratios for both observations and simulations (Figs. 3a and 8). Similar relations and seasonal variations of δ¹³C-CO₂ have been reported in other urban areas (Sturm et al., 2006; Guha and Ghosh, 2010; Moore and Jacobson, 2015; Pang et al., 2016). The simulated hourly NEE CO₂ enhancement is also shown in Fig. 7b. Note that negative values indicate net CO₂ sinks and positive values indicate net CO₂ sources. We can see large hourly variations in the growing seasons and positive enhancements during nighttime that are generally larger than negative enhancements during daytime. This shows the potential influence of NEE on δ¹³C-CO₂ seasonality. To date, no study has quantified the relative contributions to the δ¹³C-CO₂ seasonality. Here, we re-evaluate and quantify the main factors contributing to its seasonality based on the combination of δ¹³C-CO₂ observations and simulations in the following section.

Figure 7 (a) Comparisons of observed and modeled hourly δ¹³C-CO₂ from September 2013 to August 2015, where the grey line represents derived δ¹³C-CO₂ background, and (b) simulated hourly biological CO₂ enhancement. The shade and lines in both subfigures represent the periods for winter and summer, respectively.

Figure 8 Comparisons of observed and modeled (a) CO₂ mixing ratio and (b) δ¹³C-CO₂ from December 2013 to February 2014, (c) CO₂ mixing ratio and (b) δ¹³C-CO₂ from December 2014 to February 2015, (e) CO₂ mixing ratio and (f) δ¹³C-CO₂ from June to August 2014, and (g) CO₂ mixing ratio and (h) δ¹³C-CO₂ from June to August 2015.

Table 2 Statistical metrics between observed and modeled CO₂ mixing ratios and δ¹³C-CO₂ during winter, summer, and annual for 2014 and 2015. Correlation coefficient (R), averages, and root mean square error (RMSE) are displayed.

Here, we examine the comparisons for winter and summer in greater detail. The simulations showed that the model can generally capture the diel variations of observed hourly δ¹³C-CO₂ variations (Fig. 8). Statistics between observations and simulations for the two seasons are shown in Table 2. The observed seasonal average increased substantially, by 1.18 ‰, from winter 2013–2014 (−9.27 ‰) to summer 2014 (−8.09 ‰). The simulations showed a similar seasonal increase of 1.35 ‰. Some large discrepancies are evident and generally caused by the simulated total CO₂ enhancement biases (potentially caused by poorly simulated PBLH during these periods) and the negative relationship between δ¹³C-CO₂ and the CO₂ enhancement as shown in Fig. S6.

Comparisons between observations and simulations for the daily average CO₂ mixing ratio and δ¹³C-CO₂ are also shown in Fig. 9. Although the data are distributed around the 1 : 1 line for both seasons, there is less scatter and higher correlation in the winter than in the summer. We attributed this to the more complex biological CO₂ sinks in the summer, which are not adequately resolved by the relatively coarse model grid (1^∘ by 1^∘). We also performed comparisons by only choosing the daytime observations. The results indicated that daytime CO₂ mixing ratio simulations in the summer were slightly underestimated. This caused δ¹³C-CO₂ to be overestimated (Fig. S7). The simulations for winter generally captured the trends for both CO₂ and δ¹³C-CO₂ when the biological CO₂ enhancement played a relatively small role compared to anthropogenic emissions. The larger bias in the summer could result from the relatively coarse spatial–temporal resolution (aggregation error) of the Carbon Tracker biological CO₂ flux, which was 1^∘ × 1^∘ with a 3 h average. As shown in Fig. S3, the spatial distribution of land use is far more heterogeneous. This will smooth the stronger biological CO₂ signals by averaging it over the large 1^∘ × 1^∘ grid, while the urban biological CO₂ flux occurs at much finer spatial scales and likely varies at shorter time intervals.

Figure 9 Scatter plots of observed versus modeled (a) wintertime CO₂ mixing ratios, (b) wintertime δ¹³C-CO₂, (c) summertime CO₂, and (d) summertime δ¹³C-CO₂ for both years; here, these dots are daily averages.

3.2.3 Mechanisms controlling the δ¹³C-CO₂ seasonality

The mechanisms driving these seasonal variations are examined below. The peak and trough in the observed δ¹³C-CO₂ signal were observed in December and July (Fig. 10a), respectively, yielding an amplitude of 1.51 ‰. This was consistent with the simulated amplitude of 1.53 ‰. These results support the fact that the simulated δ¹³C-CO₂ seasonality agreed well with the observations (Fig. 10) and can be used to further diagnose the mechanisms contributing to the δ¹³C-CO₂ seasonality. According to Eq. (2), the δ¹³C-CO₂ seasonality can be attributed to four factors, including (1) a change in the background δ¹³C-CO₂ value from −7.64 ‰ in December to −6.66 ‰ in July, (2) a change in CO₂ background from 399 to 398 ppm, (3) the total CO₂ enhancement change from 45.7 to 37.3 ppm, and (4) the change in the isotope composition of the CO₂ enhancements causing δ_s to vary from −26.1 ‰ to −22.8 ‰.

Figure 10 Digital filtering curve fitting (CCGCRV) for background, observations, normal simulations, case 1 (excluding negative NEE when photosynthesis is stronger than respiration), and case 2 (excluding respiration and photosynthesis) in both years, (b) δ¹³C-CO₂ comparisons between normal simulations and case 1, and (c) δ¹³C-CO₂ comparisons between normal simulations and case 2.

To quantify each mechanism's contribution to the seasonality of atmospheric δ¹³C-CO₂, we recalculated δ¹³C-CO₂ by using the monthly averages as described above. First, we calculated δ¹³C-CO₂ in December and July, which were −9.54 ‰ and −8.04 ‰, respectively, with an amplitude of 1.50 ‰. Next, we replaced the δ¹³C-CO₂ background value in December (−7.64 ‰) with July (−6.67 ‰). The recalculated δ¹³C-CO₂ was −8.66 ‰ in December, indicating that the change in δ¹³C-CO₂ background value caused a change of 0.88 ‰ (9.54 ‰ minus −8.66 ‰) to the seasonality. By changing both the total CO₂ enhancement and background values, the recalculated δ¹³C-CO₂ was −8.32 ‰, contributing a 0.34 ‰ change in the seasonality (−8.66 ‰ minus −8.32 ‰). Finally, by changing δ_s from −26.1 ‰ to −22.8 ‰, together with the change in background value, the recalculated δ¹³C-CO₂ was −8.32 ‰, a change of 0.34 ‰ (i.e., −8.66 ‰ minus −8.32 ‰). This indicates that both the total CO₂ enhancement and change in δ_s contributed equally to the regional source term, causing a variation of 0.62 ‰ (i.e., 1.50 ‰ minus 0.88 ‰). Based on the above analyses, we attributed 59 % and 41 % of the δ¹³C-CO₂ seasonality to the changing δ¹³C background term and regional source terms, respectively. Further, the total CO₂ enhancement and change in δ_s, the sum of which can be treated as a regional source term, contributed equally (about 20 %) to the δ¹³C-CO₂ seasonality.

To investigate how ecosystem photosynthetic discrimination and respiration affected atmospheric δ¹³C-CO₂ seasonality, we simulated the δ¹³C-CO₂ again for two cases: (1) excluding negative NEE when photosynthesis is stronger than respiration and (2) excluding both photosynthetic discrimination and respiration. Note that only NEE was used in our study, with no partitioning between photosynthesis and respiration in the daytime. The only role of photosynthetic discrimination should be stronger than in case 1, when only negative NEE is used. The results are shown in Fig. 10b–c. Overall, the negative CO₂ enhancement caused atmospheric δ¹³C-CO₂ to become more enriched in the baseline simulations, with maximum values around 1 ‰ between April and October (Fig. 10b), and positive CO₂ enhancement (i.e., via net respiration) caused atmospheric δ¹³C-CO₂ to become more depleted compared to the baseline simulations through the whole year (Fig. 10c). By applying the CCGRCV fitting technique to the δ¹³C-CO₂ for the above two cases, we found that the δ¹³C-CO₂ seasonality decreased to 1.45 ‰ in case 1, indicating ecosystem photosynthetic discrimination explained >0.08 ‰ of the seasonality (1.53 ‰ minus 1.45 ‰). For case 2, the δ¹³C-CO₂ trough in winter slightly increased by 0.08 ‰, and the peak in summer increased by 0.20 ‰; these two factors finally led the seasonality to increase to 1.66 ‰, which was caused by much larger respiration CO₂ enhancement in summer than in winter (Fig. 7b). These results indicate that biological respiration reduced the δ¹³C-CO₂ seasonality by 0.20 ‰ and that negative NEE (photosynthetic discrimination) acted to increase the δ¹³C-CO₂ seasonality by 0.08 ‰. Generally, both ecosystem photosynthesis and respiration played minor roles in controlling the atmospheric δ¹³C-CO₂ seasonality within this urban area. In other words, the anthropogenic CO₂ emissions played a much larger role than the plants.

As shown in Fig. 5, CO₂ sources from the power industry, combustion for manufacturing, non-metallic mineral production and oil refineries and the transformation industry were the top four contributors to the CO₂ enhancements. We simulated atmospheric δ¹³C-CO₂ by assuming that no CO₂ was emitted from each of these four categories. The simulations were performed by excluding one category at a time. The results indicated that atmospheric δ¹³C-CO₂ seasonality was 1.30 ‰, 1.57 ‰, 1.30 ‰, and 1.47 ‰ when excluding the power industry, combustion for the manufacturing source, oil refineries/transformation industry, and non-metallic mineral production sources, respectively. In other words, the power industry and oil refineries/transformation industry together contributed 0.40 ‰ to the total regional source term of 0.62 ‰. The cement sources played a role in enriching by 0.07 ‰ the atmospheric δ¹³C-CO₂ in the heavy isotope, contrary to all other anthropogenic CO₂ sources.

3.3 Sensitivity analysis

3.3.1 Comparison of δ_s ⋅ ΔCO₂

Based on Eq. (2), the regional source term determines the hourly/daily variations of δ¹³C-CO₂, which is treated as a signal added to the background signal. To evaluate the model-simulated regional source term with respect to the observations, we examined daily averages for winter to minimize the influence of photosynthesis. In Fig. 11a, the observed daily δ_s ⋅ ΔCO₂ values are compared with the simulated values using the a priori anthropogenic CO₂ emissions. Here ΔCO₂ represents the total CO₂ enhancement for both observations and simulations. The product δ_s ⋅ ΔCO₂ can be interpreted as the regional source term. The average values were −1009.0 (and −841.9) ppm ‰ for observations and −1096.7 (and 1000.5) ppm ‰ for model results in 2014 (and 2015). The slope of the regression fit was 0.99 (± 0.12), and the intercept was −151.7 (± 130.1) for all data during the two winters. After applying the monthly scaling factors to constrain the anthropogenic CO₂ emissions, the re-calculated results were closer to the 1 : 1 line with a slightly improved correlation (R increased from 0.47 to 0.50; Fig. 11b). Note that the application of the monthly scaling factors only impacts the ΔCO₂ but not δ_s. The uncertainty in δ_s will be discussed next.

Figure 11 Comparisons of wintertime δ_s ⋅ ΔCO₂ using (a) a priori and (b) constrained anthropogenic CO₂ emissions.

3.3.2 Comparison between δ_{s_sim} and δ_s

To evaluate the δ_s simulations, we compared observed and simulated δ_s as displayed in Fig. 12a for all-day and nighttime conditions. Here, nighttime simulations were selected to minimize the effects of ecosystem photosynthesis and to mainly focus on the anthropogenic CO₂ sources. Two methods were used to calculate δ_s from the observations, including the Miller–Tans and Keeling plot methods. Although δ_s differed between these two methods, both displayed similar seasonal variations, with higher values (δ¹³C enrichment) in summer and lower values in winter. Such seasonal variations were also observed at other urban sites, including Beijing, China (Pang et al., 2016), Bern, Switzerland (Sturm et al., 2006), Bangalore, India (Guha and Ghosh, 2010), and Wroclaw, Poland (Górka and Lewicka-szczebak, 2013).

Figure 12 (a) Comparisons between observed and modeled δ_s, (b) relationship between cement CO₂ enhancement proportion and simulated anthropogenic δ_s for nighttime, and (c) all-day.

If the CO₂ sources/sinks are homogeneously distributed and without monthly variations, the atmospheric CO₂ enhancement components would remain unchanged, and there would be no seasonal changes in δ_s. In reality, variations in atmospheric transport processes interact with regional CO₂ sink/source changes that cause monthly variations in δ_s. The comparison of δ_s between simulations and observations indicated that the model performed well in capturing the mixing and transport of CO₂ from different sources. We can also infer from their difference that the proportions of some CO₂ categories were biased in the a priori emission map. This can be caused by both the downscaling of EDGAR inventory distribution to 0.1^∘ and the magnitude of some emissions categories. Among all anthropogenic sources, the most significant linear relations were found between the simulated anthropogenic δ_s and cement CO₂ proportions for these 24 months, with slopes of 0.33 ‰ for nighttime and 0.35 ‰ for all-day conditions (R²=0.97, p<0.001; Fig. 12b and c). These results also indicated that cement CO₂ emissions dominated monthly δ_s variations in the YRD region.

3.3.3 Sensitivity of atmospheric δ¹³C-CO₂ and δ_s to cement CO₂ emissions

The discrepancy between simulated and observed δ_s highlights that some CO₂ sources were biased in the a priori inventories. As discussed above, cement CO₂ emissions had the most distinct δ¹³C-CO₂ end-member value of 0 ‰ ± 0.30 ‰ when compared with the averages of other anthropogenic sources. Combined with its large emission compared to other regions of the world, it had a strong potential to influence δ_s and δ¹³C-CO₂. YRD represents the largest cement-producing region in the world (USGS, 2014; Cai et al., 2015; Yang et al., 2017). Its relative proportion to total national anthropogenic CO₂ emissions is about 5.5 % to 6.5 % based on the IPCC method and 7.3 % for EDGAR. These proportions are 50 % greater than the global average of 4 % (Boden et al., 2016) and much larger than most countries (Andrew, 2018) and other large urbanized areas such as California (2 %; Cui et al., 2019).

Figure 13 Sensitivity tests showing the influence of cement CO₂ emissions on δ_s for (a) nighttime, (b) all-day, and (c) the relation between cement CO₂ and δ¹³C for simulation strategies 1 (there is no bias in the total anthropogenic CO₂ enhancement such that a proportional increase/decrease in the cement component does not change the relative anthropogenic contributions) and 2 (only the cement enhancement changes). Note that the numbers in brackets indicate changes in δ¹³C with cement CO₂ enhancement proportion (the fraction of cement CO₂ enhancement to simulated total CO₂ enhancement) increase by 0.2 times. The x-axis values indicate changing cement enhancement proportions to 0.8, 1.2, 1.4, 1.6, 1.8, and 2 times the original values.

The local activity data reveal that the cement production increased from 3.55×10⁸ t in 2010 to 4.56×10⁸ t in 2014 in the YRD area. Our own calculation of the national clinker-to-cement ratio indicated a decreasing trend from 64 % in 2004 to around 56 % in 2015. Here, we applied the value of 61.7 % for 2010 and the average value of 57.0 % for 2014 to 2015. We then used the EF for clinker (0.52 ± 0.01 t CO₂ per tonne clinker; IPCC, 2019). Finally, the calculated cement CO₂ emissions were 1.14 (± 0.02) × 10⁸ t for 2010 and 1.35 (± 0.03) × 10⁸ t for 2014, indicating an 18.4 % increase over this time period. This result is close to the scaling factor of 1.145 for the total anthropogenic CO₂ emissions for the same period.

The cement CO₂ emission was 1.45 × 10⁸ t for the EDGAR products in 2010. Applying the scaling factor of 1.184, based on our independent method, the EDGAR cement CO₂ emission was 1.72 × 10⁸ t for the year of 2014. The 27 % difference between the EDGAR inventory and our independent calculations probably resulted from large errors in the clinker-to-cement ratio and regional activity data. Ke et al. (2013) reported a much higher clinker-to-cement ratio of 73 % to 70 % for China during 2005 and 2007 than the ratio of 57 % in 2014 to 2015. If we applied a 70 % ratio, the EDGAR cement CO₂ emission would change to 1.28 × 10⁸ t for 2010.

The monthly cement emission proportions varied from 6.21 % to 8.98 %, while its enhancement proportion was much larger and could reach 16.85 %. In other words, favorable atmospheric transport processes amplified the cement CO₂ enhancement proportion at our observational site (Table S2). To quantify the extent to which the cement CO₂ enhancement components can affect δ_s and atmospheric δ¹³C-CO₂, we conducted sensitivity tests by changing the cement enhancement proportions to 0.8, 1.2, 1.4, 1.6, 1.8, and 2 times its original value. These sensitivity tests are based on two different assumptions for cement CO₂ enhancement changes. (1) There is no bias in the total anthropogenic CO₂ enhancement such that a proportional increase/decrease in the cement component does not change the relative anthropogenic contributions. (2) Only the cement enhancement changes. From Eq. (2), these two assumptions will change both δ_s and δ¹³C-CO₂ but with different amplitude.

Results for the first assumption are shown in Fig. 13a–b for both nighttime and all-day δ_s simulations. The simulated δ_s increased linearly with the increase in cement proportions, at a rate of 2.73 ‰ increase per 10 % increase in cement proportions in the nighttime and 2.72 ‰ for all-day. The result for the second assumption is similar to the first one, yielding a 2.32 ‰ increase for a 10 % increase in the cement proportion. As shown in Table S2, the cement CO₂ enhancement proportions increased from 5.60 %–6.77 % (December) to 13.16 %–16.85 % (June), which is the primary cause of the observed monthly δ_s variations. The high sensitivity of δ_s to cement CO₂ proportions can partly explain the relative difference of modeled δ_s and indicates a potential advantage to constrain cement CO₂ emissions by using atmospheric δ¹³C-CO₂ observations. Finally we calculated how cement CO₂ can change atmospheric δ¹³C-CO₂ (Fig. 13c). These results show that atmospheric δ¹³C-CO₂ is more sensitive to the first assumption than the second assumption. These sensitivity analyses indicate that a cement CO₂ enhancement relative change of 20 % (or absolute 1.57 % increase) can cause a 0.013 ‰–0.038 ‰ change in the atmospheric δ¹³C-CO₂. These results indicate that δ_s is sensitive to cement CO₂ emissions.

4 Conclusions

Total annual anthropogenic CO₂ emissions for the YRD showed a high consistency between the top-down and bottom-up approaches, with a bias of less than 6 %.

Approximately 59 % and 41 % of the δ¹³C-CO₂ seasonality were attributed to the change in δ¹³C background value and the regional CO₂ source term, respectively.

The power industry and oil refinery/transformation industry together contributed 0.40 ‰ to the seasonal cycle, accounting for 64.5 % in all regional source terms (0.62 ‰).

When excluding all ecosystem respiration and photosynthetic discrimination in the YRD area, δ¹³C-CO₂ seasonality will increase from 1.53 ‰ to 1.66 ‰.

Atmospheric transport processes in summer amplified the cement CO₂ enhancement proportions in the YRD area, which dominated monthly δ_s variations. δ_s calculated from simulations was shown to have a strong linear relationship with the cement CO₂ EDGAR v4.3.2 inventory proportion in the YRD area.