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

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


Introduction
Urban landscapes account for 70 % of global CO 2 emissions and represent less than 3 % of Earth's land area (Seto et al., 2014). Such CO 2 hotspots play a dominant role in controlling the rise in atmospheric CO 2 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 2 (i.e., δ 13 C = 13 C / 12 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 2 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 2 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 2 and δ 13 C-CO 2 in urban landscapes are of great importance to monitoring these potential CO 2 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 2 emissions according to the Intergovernmental Panel on Climate Change guidelines (IPCC, 2019), and many bottom-up methods have long been used to estimate CO 2 emissions worldwide, but such methods have high uncertainties for CO 2 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 2 observations, the topdown 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 14 CO 2 -CO 2 , δ 13 C-CO 2 , and CO, is needed to better distinguish CO 2 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 δ 13 C-CO 2 observation in urban areas remains rare in inversion studies, yet such observations contain invaluable information of anthropogenic CO 2 from different categories.
Traditional estimates of δ 13 C-CO 2 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 δ 13 C-CO 2 (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 δ 13 C-CO 2 . However, most of the established long-term IRMS and IRIS δ 13 C-CO 2 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 2 and δ 13 C-CO 2 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 δ 13 C-CO 2 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 δ 13 C-CO 2 . This represents an important knowledge gap in understanding the underlying mechanisms of carbon cycling in complex urban ecosystems.
The traditional δ 13 C-CO 2 isotope partitioning methods (including the Miller-Tans and Keeling plot approaches) have been used to constrain different CO 2 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 2 enhancement components from different sources can represent CO 2 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 2 partitioning studies (Chen et al., 2017;Cui et al., 2019;Graven et al., 2018;Hu et al., 2018b). Although urban CO 2 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 2 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 δ 13 C-CO 2 and CO 2 observations. The Yangtze River delta (YRD) ranks as one of the most densely populated regions in the world and is an important anthropogenic CO 2 hotspot. Major anthropogenic sources include the power industry, oil refineries/transformation and cement production. Having the largest source of cementderived CO 2 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 2 sinks and sources within the YRD. Independent quantification of the fossil and cement CO 2 emission and assessment of their impact on atmospheric δ 13 C-CO 2 have the potential to improve our understanding of urban CO 2 cycling. Further, the observations and simulations of both atmospheric CO 2 and δ 13 C-CO 2 can help us relate atmospheric CO 2 dynamics to future emission control strategies.
Here, we combine long-term (> 2-year) CO 2 and δ 13 C-CO 2 observations with atmospheric transport model simulations to study urban atmospheric CO 2 and δ 13 C-CO 2 variations. The objectives were to (1) constrain anthropogenic CO 2 emissions and determine the main sources of uncertainty for δ 13 C-CO 2 simulations and (2) quantify the relative contributions of each factor (i.e., background, anthropogenic CO 2 emissions especially for cement production, ecosystem photosynthesis and respiration) to seasonal variations of atmospheric δ 13 C-CO 2 .
2 Materials and methods 2.1 Observations of atmospheric CO 2 mixing ratio, δ 13 C-CO 2 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 2 mixing ratios and δ 13 C-CO 2 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 2 mixing ratio and δ 13 C-CO 2 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 2 and δ 13 C-CO 2 were 0.07 ppm and 0.05 ‰, respectively. We note that the δ 13 C-CO 2 IRIS (model G1101-i) measurements are sensitive to water vapor concentration. Sensitivity tests reveal that the δ 13 C-CO 2 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). 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 2014to 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 2 -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 2 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.

10018
C. Hu et al.: Anthropogenic and natural controls on atmospheric δ 13 C-CO 2 variations 2.2 Simulation of atmospheric δ 13 C-CO 2

General equations
The simulation of atmospheric δ 13 C-CO 2 is based on mass conservation. First, we briefly describe the simulation of atmospheric CO 2 mixing ratios (more details are provided in Sect. 2.2.2), following the previous work of Hu et al. (2018b), where atmospheric CO 2 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 Note that CO 2 is the sum of all simulated sources/sinks [ CO 2_sim ] i and represents the total simulated CO 2 enhancement. We use CO 2_obs as the observed CO 2 total enhancement, which can be calculated by using the CO 2 observation minus the CO 2 background values. Based on mass conservation, we estimated the 13 CO 2 composition by multiplying the left-and right-hand sides of Eq. (1) by δ 13 C: where δ 13 C a_sim and δ 13 C bg represent the simulated atmospheric δ 13 C-CO 2 and background δ 13 CO 2 , and δ 13 i is the δ 13 C-CO 2 for end-member i (including anthropogenic and biological source categories). The δ 13 C-CO 2 contributions from all regional sources/sinks can be further reformatted as Eq. (3): 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 δ 13 C-CO 2 signals contributed by different CO 2 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 δ 13 C-CO 2 background products, and the choice of δ 13 C bg is essential for simulating δ 13 C a . Here, we apply three strategies. First, we used discrete δ 13 C-CO 2 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 δ 13 C-CO 2 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 2 background sources, most of the back trajectories originate from the free atmosphere or 1000 m higher above the ground (C. . 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 δ 13 C-CO 2 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 δ 13 C-CO 2 background when using this approach.
In the second approach, the δ 13 C-CO 2 background signal was estimated with wintertime "clean" air CO 2 and δ 13 C-CO 2 observations at the NUIST site, using the following equation: where δ 13 C a and CO 2 represent atmospheric δ 13 C-CO 2 and CO 2 observations at the NUIST site under clean conditions. Note that δ 13 C a represents the observed δ 13 C-CO 2 , not the simulated δ 13 C-CO 2 (δ 13 C a_sim ) as shown in Eq. (2).
[ CO 2_sim ] i is the simulated category-specified CO 2 enhancements. We defined clean conditions as the bottom 5 % wintertime CO 2 observations to minimize simulated CO 2 enhancement errors from both biological and anthropogenic CO 2 simulations on δ 13 C-CO 2 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 2 ] i or the corresponding δ 13 C-CO 2 signals. The hourly δ 13 C-CO 2 background value can be derived by using δ s_obs , CO 2 background, observed atmospheric δ 13 C a and CO 2 (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 2 (Alden et al., 2016;Verhulst et al., 2017), CO (Wang et al., 2010;Ruckstuhl et al., 2012) and CH 4 (Zhao et al., 2009;Verhulst et al., 2017;. To analyze the controlling factors for the δ 13 C-CO 2 seasonality, the CCGCRV (a digital filtering curve-fitting program developed by the Carbon Cycle Group, NOAA, USA) regression was applied to the back-ground, 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 δ 13 C-CO 2 .

Simulation of atmospheric CO 2 mixing ratios
In Eq. (1), the CO 2_bg is obtained from the Carbon Tracker 2016 product, which provides global CO 2 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 2 fluxes (Hu et al., 2018a, b). Considering the diurnal variations of both anthropogenic and biological CO 2 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 2 flux. The CO 2 fluxes contain anthropogenic CO 2 emissions, biological CO 2 flux and biomass burning. Here the anthropogenic CO 2 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 2 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 2 was calculated as follows: where flux i (units: mol m −2 s −1 ) corresponds to each CO 2 flux category simulated for each domain for a specific hour i, and footprint (units: ppm m 2 s µmol −1 ) is the modelsimulated sensitivity of observed CO 2 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 2 enhancements from each of the five zones were simulated by multiplying CO 2 emissions in each province by the corresponding footprint.

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 2 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. . 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 2 background of each hour, we tracked the sources of air particles' back trajectory for 7 d and defined these CO 2 mixing ratios in Carbon Tracker as the hourly CO 2 background values (Peters et al., 2007).

A priori anthropogenic CO 2 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 2 emissions ( Fig. 2a), which is based on the International Energy Agency's (IEA's) energy budget statistics and provides detailed CO 2 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 2 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 2 (ODIAC, Oda et al., 2018), also provide global CO 2 emissions. However, these inventories only provide total CO 2 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 2 category changes linearly from its 2010 value (Peters et al., 2007) and apply an annual scaling factor of 1.145 to derive CO 2 emissions for 2014 and 2015. This scaling factor is based on Carbon Tracker, dividing the same anthropogenic CO 2 emissions for the YRD in years 2014-2015 by that in 2010.
The biological flux or net ecosystem CO 2 exchange (NEE) and biomass burning CO 2 emissions come from Carbon Tracker a posteriori flux at 3 h intervals and at a spatial res- olution of 1 • × 1 • . Because NEE is much smaller than the anthropogenic CO 2 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.

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) where δ i is the δ 13 C-CO 2 value from source category i, and p i is the corresponding enhancement proportion (i.e., proportions of a specific enhancement i to total CO 2 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 δ 13 C-CO 2 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 3 plants. Since CO 2 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 3 grains that have a similar δ 13 C-CO 2 value to the biological CO 2 flux of −28.20 ‰, we assume it has the same isotope signals as local C 3 plants and ecosystem respiration. Further, the biological CO 2 flux from the Carbon Tracker assimilation system considered anthropogenic emissions to be fixed and attributed the remainder to the biological CO 2 flux (Peters et al., 2007). Consequently, we believe the uncertainty in the biological CO 2 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 13 CO 2 -CO 2 (Xu et al., 2017). We then used the results to evaluate the calculations made with Eq. (6).

Independent IPCC method for anthropogenic CO 2 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 2 emissions by using an independent IPCC method. Here, we illustrate the calculation for cement CO 2 emissions. Note that the IPCC only recommended an EF for clinker, which is an intermediate product of cement. To calculate cement CO 2 emissions, we need to calculate it based on clinker production, as shown in Eq. (7): where CO 2 [cement] is the chemical process CO 2 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 2 emission factor for clinker production. The IPCC recommended an EF clinker value of 0.52 ± 0.01 t CO 2 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.

Multiplicative scaling factor method
To quantify anthropogenic CO 2 emissions and to compare them with EDGAR products, we first derived the monthly scaling factors for anthropogenic CO 2 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 where CO 2_obs , CO 2_bio , CO 2_fire and CO 2_anthro represent observed CO 2 mixing ratios, simulated CO 2 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).

Results and discussion
3.1 Evaluation of hourly CO 2 mixing ratios

Hourly and monthly CO 2 mixing ratio comparisons
This section examines the general performance of simulating hourly CO 2 mixing ratios. The 2-year average hourly footprint is shown in Fig. 2b, where the source area (bluered) indicates strong sensitivity of the CO 2 observations to regional sources. This footprint shape is representative of the YRD area. To quantify the relative contributions from each province, we calculated CO 2 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 2 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 2 variations over the 2-year period. The model also captured the monthly and seasonal variations of CO 2 mixing ratios (daily averages are shown in Fig. S2). The simulations captured the trend of rising CO 2 mixing ratios after October and the drawdown of CO 2 to the background value during the summer. Figure 3b-d illustrate the average monthly daily, nighttime (22:00-06:00, local time), and daytime (10:00-16:00) CO 2 mixing ratios. These monthly values contain the effects of atmospheric transport, background and variations in CO 2 emissions. The observed and simulated CO 2 mixing ratios showed a significant increase from September 2013 to January 2014. Here, the CO 2 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 2 increase, and the net CO 2 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 2 mixing ratios, we compared the simulated monthly anthropogenic CO 2 enhancement differences in the same months of different years to eliminate the influence of monthly emission variations on CO 2 enhancements. Twelve monthly paired values were used and are shown in Fig. 4. This analysis indicates that atmospheric CO 2 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.
On an annual timescale, the simulated average CO 2 mixing ratios were 436.63 and 437.11 ppm for 2014 and 2015, respectively. Since the anthropogenic CO 2 emissions used in the model are the same for both years, the simulated annual average CO 2 difference can be used to quantify the influence associated with meteorological factors and ecosystem carbon cycling. Between these 2 years, the CO 2 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 2 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 2 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 2 enhancements during the non-growing seasons.

Components of urban CO 2 enhancement
Here, we diagnose the source contributions to the urban CO 2 enhancement. The observed anthropogenic CO 2 enhancements, which were derived by subtracting CO 2 background and simulated biological enhancement from CO 2 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 2 enhancements was calculated by prescribing a 2 ppm potential bias for the Carbon Tracker CO 2 fields and 50 % to the simulated biological CO 2 enhancement (Hu et al., 2018b). The corresponding simulated anthropogenic CO 2 enhancements were 36.20 and 34.90 ppm. In comparison with the simulated biological CO 2 enhancements displayed in Fig. 5a, both the observed and simulated CO 2 enhancements are indicative of a large anthropogenic (fossil fuel and cement production) CO 2 emission from the YRD.
Previous studies have also investigated urban CO 2 enhancements from a relatively broad range of developed environments worldwide. Verhulst et al. (2017) measured CO 2 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 2 enhancement was less than 20 ppm (Mueller et al., 2018). The urban CO 2 ob-  servations 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 2 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 2 enhancement, respectively. The proportions of corresponding CO 2 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 2 emissions can illustrate whether CO 2 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 2 observation site. Further, changes in source footprint caused by wind direction variations likely played an important role.

Constraints on monthly anthropogenic CO 2 emissions
To provide a robust comparison of bottom-up CO 2 emissions for the YRD, we calculated anthropogenic CO 2 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 2 emissions in 2014-2015 were 24.4 × 10 11 and 23.5 × 10 11 kg according to our own inventory and EDGAR v4.3.2 CO 2 , respectively, indicating excellent agreement (within 4 %) between these approaches. We constrained the monthly anthropogenic CO 2 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 2 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 2 emissions (Hu et al., 2018b), the monthly anthropogenic averages, and bias in negative biological CO 2 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 2 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 2 emissions were 24.6 (± 2.4) × 10 11 kg for the YRD area. The application of the MSF method provides an overall constraint on the anthropogenic CO 2 emissions (also displayed in Table 1). The main uncertainties associated with the simulation of hourly CO 2 and δ 13 C-CO 2 are uncertainty in meteorological fields, transport model (i.e., number of released particles), and a priori CO 2 fluxes. At the annual scale the main uncertainty is attributed to the PBLH simulations and a priori anthropogenic CO 2 emissions. The anthropogenic CO 2 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 2 and δ 13 C-CO 2 signals on an annual timescale.
3.2 Simulation of atmospheric δ 13 C-CO 2 3.2.1 Background atmospheric δ 13 C-CO 2 To obtain the best representative δ 13 C-CO 2 background value for the study domain, we examined the values from the three strategies described above (Fig. 6). We also compared the δ 13 C-CO 2 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 δ 13 C-CO 2 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 δ 13 C-CO 2 to be rela-tively 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 2 derived from cement production. As shown earlier, the biological CO 2 enhancement was positive in winter, which implies a positive biological CO 2 signal where ecosystem respiration is more important than photosynthesis. Further, sensitivity tests for cement CO 2 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 δ 13 C-CO 2 signal is not an ideal choice for representing the background value. The wintertime δ 13 C-CO 2 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 δ 13 C-CO 2 simulation results of Chen et al. (2006), who showed that δ 13 C-CO 2 is 0.6 ‰ higher above the PBL than in the surface layer near the ground. Recently, Ghasemifard et al. (2019) showed that hourly δ 13 C-CO 2 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 δ 13 C-CO 2 was approximately −7 ‰, during which the CO 2 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 δ 13 C-CO 2 for the study domain. Figure 7a shows the hourly δ 13 C-CO 2 simulations over a 2-year period. To the best of our knowledge, this is the first time that δ 13 C-CO 2 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 δ 13 C-CO 2 values showed seasonal variations that increased in summer and decreased in winter. This pattern mirrored the CO 2 mixing ratios for both observations and simulations (Figs. 3a and 8). Similar relations and seasonal variations of δ 13 C-CO 2 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 2 enhancement is also  shown in Fig. 7b. Note that negative values indicate net CO 2 sinks and positive values indicate net CO 2 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 δ 13 C-CO 2 seasonality. To date, no study has quantified the relative contributions to the δ 13 C-CO 2 seasonality. Here, we re-evaluate and quantify the main factors contributing to its seasonality based on the combination of δ 13 C-CO 2 observations and simulations in the following section.

Evaluation of δ 13 C-CO 2 simulations
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 δ 13 C-CO 2 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 2 enhancement biases (potentially caused by poorly simulated PBLH during these periods) and the negative relationship between δ 13 C-CO 2 and the CO 2 enhancement as shown in Fig. S6.
Comparisons between observations and simulations for the daily average CO 2 mixing ratio and δ 13 C-CO 2 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 2 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 2 mixing ratio simulations in the summer were slightly underestimated. This caused δ 13 C-CO 2 to be overestimated (Fig. S7). The simulations for winter generally captured the trends for both CO 2 and δ 13 C-CO 2 when the biological CO 2 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 2 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 2 signals by averaging it over the large 1 • × 1 • grid, while the urban biological CO 2 flux occurs at much finer spatial scales and likely varies at shorter time intervals.  Table 2. Statistical metrics between observed and modeled CO 2 mixing ratios and δ 13 C-CO 2 during winter, summer, and annual for 2014 and 2015. Correlation coefficient (R), averages, and root mean square error (RMSE) are displayed.

Mechanisms controlling the δ 13 C-CO 2 seasonality
The mechanisms driving these seasonal variations are examined below. The peak and trough in the observed δ 13 C-CO 2 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 δ 13 C-CO 2 seasonality agreed well with the observations (Fig. 10) and can be used to further diagnose the mechanisms contributing to the δ 13 C-CO 2 seasonality. According to Eq. (2), the δ 13 C-CO 2 seasonality can be attributed to four factors, including (1) a change in the background δ 13 C-CO 2 value from −7.64 ‰ in December to −6.66 ‰ in July, (2) a change in CO 2 background from 399 to 398 ppm, (3) the total CO 2 enhancement change from 45.7 to 37.3 ppm, and (4) the change in the isotope composition of the CO 2 enhancements causing δ s to vary from −26.1 ‰ to −22.8 ‰.
To quantify each mechanism's contribution to the seasonality of atmospheric δ 13 C-CO 2 , we recalculated δ 13 C-CO 2 by using the monthly averages as described above. First, we calculated δ 13 C-CO 2 in December and July, which were −9.54 ‰ and −8.04 ‰, respectively, with an amplitude of 1.50 ‰. Next, we replaced the δ 13 C-CO 2 background value in December (−7.64 ‰) with July (−6.67 ‰). The recalculated δ 13 C-CO 2 was −8.66 ‰ in December, indicating that the change in δ 13 C-CO 2 background value caused a change of 0.88 ‰ (9.54 ‰ minus −8.66 ‰) to the seasonality. By changing both the total CO 2 enhancement and background values, the recalculated δ 13 C-CO 2 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 δ 13 C-CO 2 was −8.32 ‰, a change of 0.34 ‰ (i.e., −8.66 ‰ minus −8.32 ‰). This indicates that both the total CO 2 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, Figure 9. Scatter plots of observed versus modeled (a) wintertime CO 2 mixing ratios, (b) wintertime δ 13 C-CO 2 , (c) summertime CO 2 , and (d) summertime δ 13 C-CO 2 for both years; here, these dots are daily averages.
we attributed 59 % and 41 % of the δ 13 C-CO 2 seasonality to the changing δ 13 C background term and regional source terms, respectively. Further, the total CO 2 enhancement and change in δ s , the sum of which can be treated as a regional source term, contributed equally (about 20 %) to the δ 13 C-CO 2 seasonality.
To investigate how ecosystem photosynthetic discrimination and respiration affected atmospheric δ 13 C-CO 2 seasonality, we simulated the δ 13 C-CO 2 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 2 enhancement caused atmospheric δ 13 C-CO 2 to become more enriched in the baseline simulations, with maximum values around 1 ‰ between April and October (Fig. 10b), and positive CO 2 enhancement (i.e., via net respiration) caused atmospheric δ 13 C-CO 2 to become more depleted compared to the baseline simulations through the whole year (Fig. 10c). By applying the CCGRCV fitting technique to the δ 13 C-CO 2 for the above two cases, we found that the δ 13 C-CO 2 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 δ 13 C-CO 2 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 2 enhancement in summer than in winter (Fig. 7b). These results indicate that biological respiration reduced the δ 13 C-CO 2 seasonality by 0.20 ‰ and that negative NEE (photosynthetic discrimination) acted to increase the δ 13 C-CO 2 seasonality by 0.08 ‰. Generally, both ecosystem photosynthesis and respiration played minor roles in controlling the atmospheric δ 13 C-CO 2 seasonality within this urban area. In other words, the anthropogenic CO 2 emissions played a much larger role than the plants.
As shown in Fig. 5, CO 2 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 2 enhancements. We sim- 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) δ 13 C-CO 2 comparisons between normal simulations and case 1, and (c) δ 13 C-CO 2 comparisons between normal simulations and case 2. ulated atmospheric δ 13 C-CO 2 by assuming that no CO 2 was emitted from each of these four categories. The simulations were performed by excluding one category at a time. The results indicated that atmospheric δ 13 C-CO 2 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 δ 13 C-CO 2 in the heavy isotope, contrary to all other anthropogenic CO 2 sources.

Sensitivity analysis
3.3.1 Comparison of δ s · CO 2 Based on Eq. (2), the regional source term determines the hourly/daily variations of δ 13 C-CO 2 , 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 min-imize the influence of photosynthesis. In Fig. 11a, the observed daily δ s · CO 2 values are compared with the simulated values using the a priori anthropogenic CO 2 emissions. Here CO 2 represents the total CO 2 enhancement for both observations and simulations. The product δ s · CO 2 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 2 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 2 but not δ s . The uncertainty in δ s will be discussed next.

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 Figure 11. Comparisons of wintertime δ s · CO 2 using (a) a priori and (b) constrained anthropogenic CO 2 emissions.
to minimize the effects of ecosystem photosynthesis and to mainly focus on the anthropogenic CO 2 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 (δ 13 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).
If the CO 2 sources/sinks are homogeneously distributed and without monthly variations, the atmospheric CO 2 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 2 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 2 from different sources. We can also infer from their difference that the proportions of some CO 2 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 2 proportions for these 24 months, with slopes of 0.33 ‰ for nighttime and 0.35 ‰ for all-day conditions (R 2 = 0.97, p < 0.001; Fig. 12b and c). These results also indicated that cement CO 2 emissions dominated monthly δ s variations in the YRD region.
3.3.3 Sensitivity of atmospheric δ 13 C-CO 2 and δ s to cement CO 2 emissions The discrepancy between simulated and observed δ s highlights that some CO 2 sources were biased in the a priori inventories. As discussed above, cement CO 2 emissions had the most distinct δ 13 C-CO 2 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 δ 13 C-CO 2 . 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 2 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). The local activity data reveal that the cement production increased from 3.55 × 10 8 t in 2010 to 4.56 × 10 8 t in 2014 in the YRD area. Our own calculation of the national clinkerto-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 2 per tonne clinker; IPCC, 2019). Finally, the calculated cement CO 2 emissions were 1.14 (± 0.02) × 10 8 t for 2010 and 1.35 (± 0.03) × 10 8 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 2 emissions for the same period.
The cement CO 2 emission was 1.45 × 10 8 t for the EDGAR products in 2010. Applying the scaling factor of 1.184, based on our independent method, the EDGAR cement CO 2 emission was 1.72 × 10 8 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 2 emission would change to 1.28 × 10 8 t for 2010.  . Sensitivity tests showing the influence of cement CO 2 emissions on δ s for (a) nighttime, (b) all-day, and (c) the relation between cement CO 2 and δ 13 C for simulation strategies 1 (there is no bias in the total anthropogenic CO 2 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 δ 13 C with cement CO 2 enhancement proportion (the fraction of cement CO 2 enhancement to simulated total CO 2 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 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 2 enhancement proportion at our observational site (Table S2). To quantify the extent to which the cement CO 2 enhancement components can affect δ s and atmospheric δ 13 C-CO 2 , 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 2 enhancement changes.
(1) There is no bias in the total anthropogenic CO 2 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 δ 13 C-CO 2 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 2 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 2 proportions can partly explain the relative difference of modeled δ s and indicates a potential advantage to constrain cement CO 2 emissions by using atmospheric δ 13 C-CO 2 observations. Finally we calculated how cement CO 2 can change atmospheric δ 13 C-CO 2 (Fig. 13c). These results show that atmospheric δ 13 C-CO 2 is more sensitive to the first assumption than the second assumption. These sensitivity analyses indicate that a cement CO 2 enhancement relative change of 20 % (or absolute 1.57 % increase) can cause a 0.013 ‰-0.038 ‰ change in the atmospheric δ 13 C-CO 2 . These results indicate that δ s is sensitive to cement CO 2 emissions.

Conclusions
Total annual anthropogenic CO 2 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 δ 13 C-CO 2 seasonality were attributed to the change in δ 13 C background value and the regional CO 2 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, δ 13 C-CO 2 seasonality will increase from 1.53 ‰ to 1.66 ‰.
Atmospheric transport processes in summer amplified the cement CO 2 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 2 EDGAR v4.3.2 inventory proportion in the YRD area.
Data availability. The data presented in this paper have been uploaded on our group website: https://yncenter.sites.yale.edu/ data-access (Xu and Lee, 2018).
Author contributions. CH, TJG and XL designed the study, and CH performed the model simulation and wrote the original draft. Supervision was by TJG and XL. Data acquisition was by JX, WH, DY, YC, CL, SL, and LD. All the co-authors contributed to the data analysis.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.