Partitioning of net CO2 exchanges at the city-atmosphere interface into biotic and abiotic components

Eddy covariance (EC) method has been used to measure CO2 fluxes over various ecosystems. Recently, the EC method has been also deployed in urban areas to measure CO2 fluxes. Urban carbon cycle is complex because of the additional anthropogenic processes unlike natural ecosystems but the EC method only measures the net sum of all CO2 sources and sink. This limitation of the EC method hinders us from the underlying processes of the carbon cycle, and it is necessary to partition the net CO2 fluxes into individual contributions for a better understanding of the urban carbon cycle. Here we propose a statistical method to partition CO2 fluxes into individual components by extending the method of Menzer and McFadden (2017).• Statistical method is proposed to partition CO2 fluxes into gross primary production, ecosystem respiration, anthropogenic emissions from a vehicle and building.• This method uses eddy fluxes and footprint-weighted high-resolution land cover data with temporal subsets that a few components can be negligible.• New partitioning method produces reliable individual components of the urban carbon cycle when compared to inventory data and typical biotic responses to environmental conditions.


Specifications table Subject Area
Earth and Planetary Sciences More specific subject area Urban Climate, Microclimate, Flux measurement, Carbon Cycle Method name CO 2 flux partitioning in urban area Name and reference of original method Crawford and Christen [1] . Spatial source attribution of measured urban eddy covariance CO 2 fluxes. Theoretical and applied climatology, 119 (3)(4), 733-755. Menzer and McFadden [12] . Statistical partitioning of a three-year time series of direct urban net CO2 flux measurements into biogenic and anthropogenic components. Atmospheric Environment, 170, 319-333. Resource availability * Method details Urban CO 2 budget is formulated as: (1) where F C is the net CO 2 exchange at the city-atmosphere interface, RE is the ecosystem respiration, P is the gross primary production by vegetation, and E R and E B are the anthropogenic CO 2 emissions from fossil fuel combustions by vehicles and buildings, respectively. It is important to partition F C into individual components (i.e., RE, P, E R , and E B ) to investigate the controlling factors of the urban carbon cycle in a changing environment. In natural canopies and croplands, RE is estimated by modeling the relationship between the EC-measured nighttime F C and air or soil temperature, and then P is estimated by subtracting this modeled RE from F C (e.g., [15] ).
However, unlike natural ecosystems, there are extra anthropogenic CO 2 sources (i.e., E R and E B ) at the city-atmosphere interface. Thus, additional information and process are required for the F C partitioning into each CO 2 sources and sink. For example, Pataki et al. [13] used stable isotope measurement data for the partitioning. Meanwhile, Crawford and Christen [1] and Menzer and McFadden [12] proposed statistical methods for the F C partitioning in urban areas. Crawford and Christen [1] partitioned F C by using meteorological variables and the flux footprint weighted land cover climatology in suburban area where vehicle emission is a dominant source. Menzer and McFadden [12] partitioned F C by modeling the anthropogenic components with traffic volume and air temperature ( T air ) based on winter F C data when only anthropogenic sources contribute to net CO 2 fluxes. Our new proposition extends the statistical method of Menzer and McFadden [12] (MM2017 hereafter) by combining MM2017 with the flux footprint weighted land cover climatology based on Crawford and Christen [1] for the urban CO 2 flux partitioning. Similarly to MM2017, our method selects specific temporal subsets in that some components in Eq. (1) are negligible. Additionally, our method takes advantage of traffic volume data and the footprint weighted road cover fraction ( λ R ). Accordingly, our proposed method is applicable to more general urban conditions especially when E R is the main source from the heavy traffic volume in a city. Fig. 1 shows the schematic flow chart to partition net CO 2 exchange into contributing components by our method with the MM2017 for comparison. Here we explain the procedures of our method with an example EC data observed in an artificially constructed urban forest [4 , 8 , 9] .

Data preparation
Our partitioning method requires the EC data with other ancillaries of the footprint-weighted land cover fraction ( λ), traffic volumes, and vegetation index (e.g., EVI: Enhanced Vegetation Index). λ is estimated by overlaying flux footprints on the high-resolution land cover map, which needs the atmospheric stability estimate from the EC data, a footprint model, aerodynamic roughness parameters, and the morphological data of digital surface and terrain models.
In the example case, the 30-min CO 2 flux is computed using the EddyPro software (version 6.2.0, LI-COR, USA) from 10 Hz sampled raw data and the post data process for quality control ( [3] and references therein). Footprint climatology is computed using the footprint model of Hsieh et al. [2] and a 1 m-resolution land cover map provided by Environmental Geographic Information  [12] to partition CO 2 flux ( F C ) into four components: ecosystem respiration ( RE ), gross primary production ( P ), and anthropogenic CO 2 emissions by vehicles ( E R ) and buildings ( E B ). The hatted variables are the modeled components in Eq. (1) . The data filter conditions used in the process of modeling each component are not universal, and it is recommended to set the criteria according to the characteristics of each observation site. Table 1 Comparison between method of (a) the proposed method in this study and (b) Menzer and McFadden [12] . The parameters used to model each component are values that are only applicable to their corresponding site.

Partitioning method (a) The proposed method (b) Menzer and McFadden [12]
Partitioning order Service ( https://egis.me.go.kr ) following Kim et al. [6] . In the land cover map, the land cover was classified into tree, grass, building, road, impervious, water, bare soil, and the other. As inputs for the footprint model, aerodynamic roughness parameters such as roughness length and zero-plane displacement height are estimated using the morphological data of digital surface and terrain models provided by National Geographic Information Institute ( https://www.ngii.go.kr ) [4 , 8] at 1 × 1 m 2 resolution. Eventually, λ is computed by the summation of each footprint-weighted land cover [1] .
To estimate E R over the heterogeneous urban area, E R is modeled with traffic volumes and footprintweighted road cover fraction ( λ R ) in our method. This procedure presumes that E R is proportional to its corresponding surface type fraction (i.e., λ R ). In the example, the hourly traffic volume is counted on a road near the flux tower ( https://topis.seoul.go.kr ).

Partitioning procedures
In statistical partitioning methods, it is important to select relevant temporal subsets of F C data so that components in Eq. (1) are negligible. For example, MM2017 chose wintertime F C with low traffic volume in cold temperate climate zone so that P and RE were negligible and E R was also small, thus only E B contributed to F C (Step 1 in Fig. 1 a). However, in a typical metropolis, such condition is not always feasible because it is probable that traffic volume is high in the footprint areas. Our method is designed to estimate E R in such complex conditions ( Table 1 ). The following sections describe the procedure to partition RE, E R , E B , and P in detail.

Estimation of ecosystem respiration
First, RE is computed by the van't Hoff-type regression of nighttime F C to T air [17] because RE has the exponential relationship with T air generally. This nonlinear regression should be limited to the periods of negligible P, E R , and E B : 1) P is zero during nighttime, 2) E R is negligible if the source area weighted road fraction is small (i.e., λ R ∼ = 0), and 3) E B is insignificant in the warm season because people stop the use of gas-based heating systems.
In the example case, nighttime data is selected if downward shortwave radiation ( K ↓ ) is below 5 W m −2 during the sunset periods, and the warm season is defined as the period of T air > 12.0 °C which is determined by the changing point analysis from the exponential decreases to increases of F C  to T air based on Killick et al. [5] ( Fig. 2 ). The final estimation of RE is as follows ( Fig. 3 ): Q 10 (the rate by which respiration is multiplied when temperature increases by 10 °C) is approximately 1.9 at the example site and matched the Q 10 value for ecosystem respiration (2.2 ± 0.7) calculated in natural forests across 42 FLUXNET sites [10] .

Estimation of anthropogenic CO 2 emissions from fossil fuel combustions by vehicles
E R is estimated by the data during the transition period (hereafter TP) defined by the period of both negligible E B and P . TP corresponds to the leaf-off period with the relatively warm temperature and no space heating. That is, during TP, F C consists of only RE and E R and so E R is simply obtained by subtracting the modeled RE (i.e., Eq. (2 )) from F C .
In the example case, the leaf-off period (1st November -14th April) is decided based on EVI, and TP is selected as total 25 days for two years with daily mean T air > 12.0 °C in this leaf-off period. In previous studies, only traffic volume was used as an input variable to predict E R , but E R shows diurnal variation concurrently with traffic volume and λ R ( Fig. 4 a-c). Indeed, E R shows the best correlation with traffic volume multiplied by λ R . We binned the data by averaging every 20 samples in ascending order from the lowest traffic volumes multiplied by λ R and regressed them against E R with the atmospheric stability ( Fig. 4 d):

Estimation of anthropogenic CO 2 emissions from fossil fuel combustions in buildings
After estimation of RE and E R , E B can be estimated by subtracting them from the observed F C during the leaf-off season (i.e., P ∼ = 0 ). E B consists of emissions from space heating ( E B(s) ) and from water heating ( E B(w) ). In general, E B(s) is proportional to T air because natural gas consumption for space heating increases as it gets colder, and E B(w) is small in nighttime. Accordingly, nighttime E B data consists of E B(s) only and E B(s) is estimated from the relationship between nighttime E B and T air . In the example case, nighttime E B(w) is also small because the urban park is closed at night (from 21:00 to 08:00), thus leading to no usages of hot water. Particularly, in the example case, E B(s) is estimated separately for the first and second year, and on the western and the eastern side, respectively to incorporate spatio-temporal variations in E B(s) around the measurement location ( Fig. 5 ): ( eastern side , 2 nd year ; p < 0 . 0 0 01 ) (9) The slopes in Eq. (6) -(9) are the sensitivity of E B(s) to T air and their annual changes match variations in annual gas consumption recorded at the park facility building around the tower. In the example, a direct emission source by space heating is located on the western side of the tower and is larger in the second year than the first year based on the records at the facility building around the tower. Notably, these annual changes in the slopes match the larger slope in Eq. (7) compared to that in Eq. (6) ( Fig. 6 ).
There is still non-zero residual during the leaf-off season with its diurnal variation, which can be attributed to E B(s) and errors. In our example, on the western side, it shows a clear diurnal pattern and follows the number of park visitors but the eastern side residuals change their sign for a day ( R 2 = 0.74, Fig. 7 ). Our results suggest that this residual mainly comes from E B(s) in daytime but error is dominant in case of no significant CO 2 sources (i.e., the eastern side in the example). Furthermore, the monthly sum of the modeled E B is comparable to the estimation based on the monthly gas consumption data of the park facility during the winter with a significant correlation ( R 2 = 0.94), indicating that our method is more reliable in monthly timescale ( Fig. 6 ).

Estimation of gross primary production
After RE, E R and E B are modeled sequentially, P is obtained as the residuals by subtracting the remaining components from the observed F C ( Fig. 8 ):  . In this step, it is prerequisite to fill the missing in F C and small negative P ( < ∼ 2 μmol m −2 s −1 ) at night is replace with zero. Further investigation reveals that 10% errors in turbulent fluxes and their related uncertainties in the flux footprint produces P uncertainty of < 1% only.

Concluding remarks
This study proposes a new statistical methodology for separating the eddy covariance CO 2 flux observed in a complex urban environment with air temperature, traffic volumes, and land cover information. Our proposed method can be used for urban interfaces influenced by heterogeneously distributed CO 2 sources and sink where previous statistical methods are inapplicable. Our method is applied for an example data observed at an urban forest where roads are adjacent to the measurement location and provides comparable estimation to the inventory data. In typical urban environment where CO 2 sources and sinks have spatiotemporally complex distribution, our methodology can provide reliable estimations of CO 2 sources and sink and improve our understanding of urban carbon cycle accordingly and our estimations are insensitive to uncertainties in measurement and ancillary data. More extensive observation and validation studies are required to investigate impacts of uncertainties in the footprint model and EC data on the flux partitioning. It is also notable that our estimated empirical parameters are site-specific and need independent estimation in other urban sites.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.