Investigating the Uncertainties Propagation Analysis of CO₂ Emissions Gridded Maps at the Urban Scale: A Case Study of Jinjiang City, China

Abstract

Gridded CO₂ emission maps at the urban scale can aid the design of low-carbon development strategies. However, the large uncertainties associated with such maps increase policy-related risks. Therefore, an investigation of the uncertainties in gridded maps at the urban scale is essential. This study proposed an analytic workflow to assess uncertainty propagation during the gridding process. Gridded CO₂ emission maps were produced using two resolutions of geospatial datasets (e.g., remote sensing satellite-derived products) for Jinjiang City, China, and a workflow was applied to analyze uncertainties. The workflow involved four submodules that can be used to evaluate the uncertainties of CO₂ emissions in gridded maps, caused by the gridded model and input. Fine-resolution (30 m) maps have a larger spatial variation in CO₂ emissions, which gives the fine-resolution maps a higher degree of uncertainty propagation. Furthermore, the uncertainties of gridded CO₂ emission maps, caused by inserting a random error into spatial proxies, were found to decrease after the gridding process. This can be explained by the “compensation of error” phenomenon, which may be attributed to the cancellation of the overestimated and underestimated values among the different sectors at the same grid. This indicates a nonlinear change between the sum of the uncertainties for different sectors and the actual uncertainties in the gridded maps. In conclusion, the present workflow determined uncertainties were caused by the gridded model and input. These results may aid decision-makers in establishing emission reduction targets, and in developing both low-carbon cities and community policies.

1. Introduction

Cities play an essential role in low-carbon development [1] since they accommodate half of the global population, and account for three quarters of both global energy consumption and greenhouse gas emissions [2]. Understanding the spatiotemporal distributions of CO₂ emissions from anthropogenic activities is essential for their mitigation [3,4]. Geospatial datasets obtained from multiple sources (e.g., remote sensing satellites, crowd-sourcing data, location-based service data, and sensors data) are increasing in number, which enables the accurate measurement of the spatiotemporal dynamics of CO₂ emissions. However, large uncertainties in estimating CO₂ emissions at the urban scale, reduce the efficacy of mitigation policies [5]. For instance, uncertainties (50–100%) at the state and urban scales in the USA [6] are greater than the emission reduction target for New York City [7]. High uncertainty (42.96%) has also been reported between consumption-based (142 Mt) and production-based emissions (81 Mt) in Beijing, China [8]; with a value nearly equal to the emission reduction target [9]. Uncertainties for medium and small cities in China are generally greater than 40 and 20%, respectively [8].

Estimating CO₂ emissions at the urban scale, including both un-gridded (i.e., administrative unit maps[10,11,12,13,14]) and gridded maps, cannot avoid the propagation of uncertainties from input to result, which highlights the importance of being aware of uncertainty estimation, especially in gridded maps due to its implications for the precision mitigation of CO₂ emissions. Two primary sources of uncertainties for the CO₂ emission gridding process can be identified as the uncertainties caused by the gridded model (geoprocessing analysis), and those caused by the input. Previous studies have focused on the uncertainties caused by the input. For example, Gately et al. (2017) investigated the two detailed main sources of uncertainties for input: (1) the total emissions estimations, and (2) the spatial (or temporal) distributions of proxies (i.e., proxies refer to the potential geospatial datasets that could reflect and be used to generate the spatial or temporal distribution of CO₂ emissions) [5]. Uncertainties in total emission estimations have been analyzed and reported in several studies on un-gridded maps [12,15]. These uncertainties have been attributed to activity levels (energy consumption or industrial production, e.g., statistical data of energy consumption) and emission factors [16]. A previous study reported a 20% difference between the aggregated statistical data on energy in 30 provinces, and the national statistical data on energy for China [16]. Moreover, the default global emission factors do not account for the actual emission scenarios among different countries, and tend to overestimate the total emissions when compared to local reports (e.g., a 14% emission bias reported in China) [17]. Furthermore, uncertainties in the total emission estimations vary at different scales. Previous studies have reported low uncertainty (−6 to 6%) in the total emissions estimations for 24 Chinese cities in 2010, medium uncertainty (−9 to 11%) in the national CO₂ emissions estimate for China in 2005, and high uncertainty (−18.75 to 18.75%) in the global CO₂ emissions estimate from combustion sources in 2007 [15,18,19]. For the uncertainties in the spatial (or temporal) distributions of proxies, these uncertainties mostly occur in gridded maps, and can be attributed to the aggregated effects of their spatial resolutions, and the inaccurate spatial patterns of emissions [5]. The aggregated effects misestimate emissions at the grid-cell level, because emissions are generally generated at a spatial scale coarser than that of the actual source activities. For example, a comparison of the Fossil Fuel Data Assimilation System, the Open-Source Data Inventory for Anthropogenic CO₂, the Emissions Database for Global Atmospheric Research (EDGAR), and the Anthropogenic Carbon Emissions System indicated significant differences in CO₂ emissions at regional (20%) and city scales (50–250%) [5]. Additionally, the inaccurate spatial patterns of emissions of spatial proxies mean that the uncertainty exists in a digital representation (e.g., uncertainties in retrieval products, uncertainties in geoprocessing process, or the low-quality satellite imagery due to cloud mask), or that spatial proxies cannot entirely reflect the spatial distribution of CO₂ emissions (e.g., population density gridded maps) [20]. For instance, the gridded uncertainties in the CDIAC data were based on uncertainties in the spatial proxies, which were attributed to the transformation of the coordinate system (29–112 km per degree at different latitudes, 0–73.87%) and the raster representation (0–100%). The results indicated that uncertainties in the CDIAC map ranged from 4 to 190% [21].

However, to our knowledge, few studies have investigated the uncertainties caused by the gridded model. Most of the existing detailed gridded CO₂ emission maps have been classified into two categories, namely the “global downscaled” and “bottom-up” approaches [22]. The global downscaled approach uses spatial proxies, such as population data and nighttime light imagery, to determine the total global, national, and regional emissions over a defined space and time domain [23]. Whereas, the bottom-up approach focuses on spatiotemporally explicit gridded CO₂ emission maps based on direct flux monitoring (eddy covariance measurements) and sector emission data from various sources, such as transportation, industry, and commercial buildings [22]. In both the allocated model of the global downscaled approach, and the transformed model of the bottom-up approach, may exist model structural error, leading to uncertainties. Furthermore, uncertainties existing in non-spatial and spatial data will subsequently propagate through the gridding process of CO₂ emissions. Generally, the uncertainties propagated through the gridding process are ultimately enlarged/reduced. This may cause misleading findings for CO₂ emissions-related studies.

Thus, the quantification of uncertainties at the urban scale is central to research on CO₂ emissions [5]. Previous methodologies for the quantification of uncertainties of CO₂ emissions could be categorized into two main types: using Monte Carlo simulations, or using the inventory comparison method. Monte Carlo simulations can be used to estimate uncertainties in total emission estimations by generating samples according to the probability distribution function, and prior coefficients of variance for activity levels and emission factors [15,18,19]. Previous studies reported prior coefficients of variance for activity levels, ranging from 5 to 30% for different sectors [15,19]. The coefficients of variance for emission factors were primarily based on Intergovernmental Panel on Climate Change guidelines [24]. The existing Monte Carlo simulations can estimate the total uncertainties but lack the spatial distribution of uncertainties to enable identifying both the high-emission and low-uncertainty areas. The inventory comparison method is a commonly used method for evaluating the uncertainties in spatial (or temporal) distribution of proxies by comparing the gridded maps with other open products. The inventory comparison method [5] can be used to evaluate uncertainties in the sum of all cell values, and to subsequently draw uncertainty distribution maps. There is a study that compared the National Bureau of Statistics of China, EDGAR, and Fossil Fuel Data Assimilation System with the China High Resolution Emission Database (CHRED) at two levels, and demonstrated strong correlations (correlation coefficients of 0.58–0.86) between CHRED and the other products [25]. The total emissions estimate obtained from high-resolution mapping methods (7.83 Pg·C·yr⁻¹) was similar to that obtained by the International Energy Agency (7.873 Pg·C·yr⁻¹) [18]. The inventory comparison method has limited applications because of the lack of emission inventories at the urban scale. Furthermore, different accounting scopes among the open gridded CO₂ emission inventories can also cause incorrect uncertainty evaluation results.

To address the above mentioned issue, an analytic workflow was proposed to analyze the propagated uncertainties caused by the gridded model and the input for gridded CO₂ emission maps. The present workflow used four sub-modules based on Monte Carlo simulations and a bootstrap sampling method to analyze uncertainties, without other detailed open emission inventories. Two of the submodules obtained the corresponding uncertainty of each grid value, generated the uncertainty map caused by the gridded model, and the uncertainties for the sum of each cell (also referred to as propagated uncertainties caused by the model); the other submodules generated the uncertainty distribution maps based on the total emission estimations and spatial proxies (also referred to as propagated uncertainties caused by input). We regarded the gridded maps of CO₂ emissions constructed in previous studies as a case study, and applied the workflow to estimate the uncertainties [26,27,28,29]. The estimation of different uncertainties helps decision makers in formulating relevant policies. Uncertainties in total emission estimations aid the determination of emission reduction targets, the corresponding risks for cities and enterprises, and significant emission sources. Uncertainty maps could help to identify locations suitable for developing low-carbon communities. Therefore, this workflow can be adopted in other studies which focus on urban-scale CO₂ emissions. Additionally, the estimation of different sources of uncertainty could improve emission-reduction policies.

2. Materials and Methods

2.1. Overview

Taking Jinjiang in 2013 as the case study, gridded maps of annual CO₂ emissions were produced with sizes of 30 m and 500 m using two multi-sources geospatial datasets of corresponding resolution, including satellite (e.g., night light imagery), vector data (e.g., road network, town-level populations), and satellite-derived products (e.g., digital elevation model). Our goal was to quantify the uncertainties caused by the gridded model, and inputs of the gridded CO₂ emission maps at the urban scale. We proposed an analytic workflow based on the Monte Carlo simulation and bootstrap sampling method. The workflow included four different submodules, described in Section 2.3.

2.2. The Description of the Case Data

The CO₂ emission gridded maps of Jinjiang (E 118°24′–118°43′, N 24°30′–24°54′), a southeastern city in China (Area 649 km²), were constructed based on geospatial datasets (e.g., remote sensing satellite-derived products) from previous studies [26,28]. Jinjiang is one of the most developed counties in China; yet, it represents the local urbanization mode, with fast economic growth and urban sprawl, a stable permanent population, and a large floating population. The gridded maps were constructed using a general hybrid approach, based on global downscaled and bottom-up elements (e.g., industrial areas). In the previous studies, the first step calculated the total (sum of emissions from three sectors below), residential, industrial (except the energy consumption of the electrical production department), and transport emissions using energy consumption values within the urban geographic boundary of Jinjiang City in 2013, and the corresponding emission factors. The next step was generating three corresponding spatial proxies of each sector: resident (refers to population density gridded maps), industry (refers to the product of a binary-layer of industrial land and night light intensity), and transport (refers to the road area). These spatial proxies were produced at 30 m and 500 m, and the sum of all the grid values was kept as 1. Each grid value, usually named “weight”, represents the ratio of the corresponding grid occupied with the specific emissions. Finally, using the following formula the total gridded CO₂ emission maps were generated. The coarser-resolution gridded CO₂ emission maps were not aggregates of the fine-resolution gridded CO₂ emission maps. They were produced by using two-resolution geospatial datasets.

$$Gri{d}_i={\displaystyle \sum }_{l=1}^3{C}_l·Weigh{t}_{i,l}={\displaystyle \sum }_{l=1}^{3}A{L}_l·E{F}_l·Weigh{t}_{i,l}$$

(1)

where Grid_i is the total CO₂ emissions (unit: t) at the ith grid (i = 1, 2, 3..., n), C_l (units: t) is the total amount of CO₂ emissions from different sectors, AL_l (units: t) is the total energy consumption from different sectors, and EF_l (units: t/t CO₂) is the emission factor for different sectors based on the IPCC method [24] at the ith grid (l = 1, 2, 3, which represent residential, industrial, and transport emissions, respectively), and Weight_i,l is the weight of the specific sector on the grid, i. In fact, Weight_i,l is the mathematical form of the spatial proxies. Detailed information about these gridded maps is shown in the paper of Dai et al. (2020) [26].

2.3. Analytic Workflow of Uncertainty Propagation

An analytic workflow including four submodules was developed to quantify the uncertainties caused by the gridded model and input in the gridded maps of CO₂ emissions (Figure 1). The uncertainty results of submodules were U1, U2, U3, and U4.

2.3.1. The Design of Analytic Workflow

The initial step in estimating the uncertainties propagated during the gridding process was analyzing the gridded model and inputs, to understand the potential uncertainty sources and the propagating process. In Figure 1, the gridding process consists of raw data (e.g., energy, industrial points, roads, and night light imagery satellite), standard input data (e.g., activity levels, emission factors, and spatial proxies), and the gridded model (i.e., a raster-based model to overlay the gridded CO₂ emission maps of three sectors using Formula (1)). Both the raw data and standard input data belong to the input of the gridded model. Hence, the potential uncertainty sources, which may be propagated through the gridding process, include the uncertainties caused by inputs, namely uncertainty 1, and the uncertainties caused by the gridded model, namely uncertainty 2. The gridded model causes the propagation of uncertainties 1 and 2. Overall, uncertainties of output depend on the synergistic effect of uncertainty 1, uncertainty 2, and uncertainty propagation.

The four submodules of the workflow were based on the two basic assumptions: (1) the input is “true” and has not caused the uncertainties, and (2) the gridded model has not caused the uncertainties.

• Assumption 1: the input is “true” and has not caused the uncertainties

We developed the first two submodules (Submodules 1 and 2) based on assumption 1. In these submodules, the uncertainty of output only depended on the uncertainties caused by the gridded model and uncertainty propagation.

For the first submodule, the global probability density functions (PDFs) of “true” spatial proxies among three sectors were fitted using a maximum likelihood estimate method, and we verified the goodness of fit using a Kolmogorov–Smirnov test (KS test) [30]. Once the KS test was significant, the fitted PDFs would be used to generate simulated samples using Monte Carlo simulation. Subsequently, we compared the cumulative distribution functions (CDFs) of the simulated and “true” grid values. The differences between the two CDFs represented the grid value–uncertainty relationship in the final gridded CO₂ emissions maps caused by the gridded model. Then the look up table (LUT), a table including each grid value and its corresponding uncertainty, was generated by the above CDFs comparisons. Finally, the uncertainty gridded maps caused by the gridded model were generated by the corresponding relationship between grid value and uncertainty, from the LUT and gridded CO₂ emission maps. These uncertainties were named U1.

For the second submodule, a bootstrap sampling method [31] was used to generate the uncertainty for the sum of the grid values caused by the gridded model. Then, we randomly selected 1000 samples from the grid values of three spatial proxies among three sectors, and repeated the sample generation process (sample with replacement) 10,000 times. Once the sample distributions were obtained, we calculated the mean grid value of total, residential, industrial, and transport CO₂ emissions for each sample generation process. Using the samples of mean grid values, we could estimate the mean value, and the 95% confidence interval (CI), of the sum of the grid values among total, residential, industrial, and transport gridded CO₂ emission maps through multiplying the number of grids; this was named U2.

• Assumption 2: the gridded model has not caused the uncertainties

We developed other submodules (Submodule 3 and 4) based on assumption 2. In these submodules, the uncertainty of output only depended on the uncertainties caused by the input and uncertainty propagation. The input could be categorized into two detailed classifications: total emissions estimation, and spatial proxies. Submodules 3 and 4 aimed to investigate the above corresponding uncertainties caused by the total emission estimation and spatial proxies, respectively.

For the third submodule, Monte Carlo simulation was used to estimate the uncertainty maps caused by the total emission estimations (i.e., activity levels and emission factors). In this study, we took a city, just as an experiment, to investigate the uncertainties caused by the total emission estimations. Hence, we only considered the uncertainty of activity levels. Regarding the total CO₂ emissions, we assumed that the emission factors and spatial proxies were accurate, and that the uncertainties in activity levels followed a normal distribution [18,19]. Once the coefficient of variation of the activity levels had been obtained, Monte Carlo simulation was used to generate 10,000 samples of activity levels in order to generate new gridded maps of CO₂ emissions. U3 was obtained as the difference between the new gridded maps and the “true” gridded maps (our gridded maps generated in the previous studies). Furthermore, submodule 3 could also be used to estimate the uncertainties caused by emission factors.

For the fourth submodule, a random error of the spatial maps was used to produce the uncertainty maps caused by the spatial proxies. Specifically, we randomly selected 20% of the grids, and set each of these grid values as a sum of the mean value of the selected grids and the random error, which was generated by a normal distribution with a zero mean and fixed standard deviation (i.e., the mean value of the selected grids divided by six). Then we used the new spatial proxies with random errors to generate the new gridded maps of CO₂ emissions. We compared these new gridded maps with the “true” gridded maps, and thus obtained U4.

It should be noted that no field measurement data were used in this study, and the uncertainties of the gridded CO₂ emissions maps were the expected results based on probability theory. Detailed formula and experiment settings will be introduced in Section 2.3.2.

2.3.2. The Formula of the Submodules and Experiment Setting

All the grid values of the CO₂ emission maps were the stochastic variables, Grid_i, within a product of the total CO₂ emissions estimate, C_l, and the spatial proxies, Weight_i,l, of an unknown probability distribution, which satisfy Equation (2).

$$\left\{\begin{array}{c}Gri{d}_i={{\displaystyle \sum }}_{l=1}^3{C}_{l}Weigh{t}_{i,l}={{\displaystyle \sum }}_{l=1}^{3}A{L}_l·E{F}_l·Weigh{t}_{i,l}\\ Weigh{t}_{i,l}~P_{un,l}\left({\lambda }_1,\cdots ,{\lambda }_i\right)\end{array}\right.$$

(2)

where the P_un,l(…) was the lth probability density function of the spatial proxies. The spatial proxies obeyed the unknown probability distributions which could be fitted. λ₁, …, λ_i were the fitted parameter of probability distributions of the spatial proxies. The descriptions of Grid_i, C_l, AL_l, EF_l, and Weight_i,l are the same as the descriptions of Equation (1). We developed and conducted the four submodules based on the above theory. Detailed formulas and experimental settings of the submodules are mentioned below.

Submodule 1 indicated the difference between the CDFs of the simulated and actual grid values, which elucidated the relationship between each cell value and its corresponding uncertainty. First, the global PDFs of the spatial proxies for CO₂ emissions among three sectors at 30 and 500 m were fitted by the maximum likelihood estimation method, and the KS test was performed to verify the fitting results [30]. As the PDFs of the spatial proxies were similar to the logarithmic normal distribution, we assumed that the PDFs of the spatial proxies were this distribution, and fitted the PDFs. Six global PDFs from three sectors at two resolutions for the entire area were fitted. Based on six fitted PDFs, Monte Carlo simulation was used to generate 100,000 samples for each sector at two resolutions. Using Formula (1), we could obtain the total, residential, industrial, and transport CO₂ emissions at each grid value. Hence, we obtained the CDFs of the simulated grid values and actual values of CO₂ emissions, which were then used to compare the simulated and actual grid values. Defining the uncertainty as the difference among two CDFs, and two LUTs (30 and 500 m), the grid value and its corresponding uncertainty were generated by self-defined intervals of grid values. Finally, we retrieved and generated the uncertainty maps from the gridded model by searching the corresponding uncertainty of each grid value in LUT.

Submodule 2 evaluated uncertainty in the sum of the grid values (i.e., total emission estimations) by the non-parametric bootstrap sampling method [31]. During each sampling process, a total of 1000 grid values of CO₂ emissions, based on each of the three spatial proxies at 30 and 500 m, were sampled, and the sampling process was repeated 10,000 times. Subsequently, the 1000 corresponding grid values of CO₂ emissions were calculated using Formula (1). Thus, we obtained the mean values of 1000 corresponding grids in each sampling process. Then 10,000 mean values of CO₂ emissions ($\overline{Gri{d}_1},\overline{Gri{d}_2},\cdots \overline{Gri{d}_m}$ m = 1, 2, …, 10,000) in the sampling grids were obtained, which could represent the probability distribution of mean values for all the grids ($\overline{Gri{d}_i})$, due to the law of large numbers. Hence, the PDFs, statistical terms, and the 95% CI of total CO₂ emissions for all the grids were estimated through linear transformation. Since total CO₂ emissions for all the grids were equal to the product of the mean value and the number of grids, the 30 m and 500 m maps had 1,457,372 and 5332 grids, respectively.

Submodule 3 was used to evaluate the propagated uncertainty caused by the total emission estimations. It was based on the generative uncertainty of the total emission estimations, which was in turn based on the coefficients of variance of emission factors, activity levels, or both. These were used to generate normally distributed samples, and to calculate their means and CI. Zhao et al. [19] reported a 5–30% variation in the normal distributions of activity levels for different sectors of CO₂ emissions. In this study, the coefficients of variance for activity levels of residential, industrial, and transport emissions were set to 20, 10, and 16%, respectively [19]. Monte Carlo simulations were used to generate 10,000 activity level samples for each sector, which generated 10,000 new gridded CO₂ emissions maps. The percentiles for CO₂ emissions at each grid were obtained from the new gridded maps. Therefore, uncertainty and percentile values were calculated as follows:

$$\left\{\begin{array}{c}AL~N\left({\mu }_{AL,l}, {\sigma }_{AL,l}{}^2\right), cv=\frac{{\sigma }_{AL,l}{}^2}{{\mu }_{AL,l}}\\ Uncertaint{y}_i=\frac{C{I}_{i,95}}{{\overline{CO}}_{2,i,N}}\\ C{I}_{i,95}=\frac{C{O}_{2,\mathrm{i},97.5}-C{O}_{2,\mathrm{i},2.5}}{2}\end{array}\right.$$

(3)

where AL obeys the normal distribution with μ_AL,l and σ_AL,l², cv is the assumed coefficients of variance for activity levels and equal to the σ_AL,l² divided by μ_AL,l, μ_AL,l is equal to the total CO₂ emissions of the lth sector, σ_AL,l² is calculated by μ_AL,l and cv, Uncertainty_i is the uncertainty in CO₂ emissions at the ith grid,${\overline{CO}}_{2,i,N}$is the mean CO₂ emissions at the ith grid obtained from Monte Carlo simulations, CI_i,95 is the 95% CIs width of CO₂ emissions at the ith grid in the Monte Carlo simulations, and CO_2,i,97.5 and CO_2,i,2.5 are the 97.5 and 2.5% quantiles of CO₂ emissions at the ith grid in the Monte Carlo simulations, respectively.

Submodule 4 evaluated the propagated uncertainty caused by spatial proxies, through inserting a random error into the spatial proxies or quality assessment. Most spatial proxies are difficult to evaluate the accuracy of without validation datasets. Hence, the prior error of the spatial proxies will be used to investigate the propagated uncertainty caused by spatial proxies. Furthermore, some spatial proxies caused the evaluated uncertainties during the generation process, such as the quality assessment files of the satellite-derived products. In this study, 20% of the grids of the spatial proxies were randomly sampled. These grid values were re-generated using a sum of the mean value of the sampled grids and a random error, ξ_i. The random error obeys the normal distribution, whose mean is equal to 0 and standard deviation is equal to the square of the mean value of the sample grids divided by 3. The setting of the random error considered two limited rules: (1) the sum of the spatial proxies must be 1, and (2) each grid value of the spatial proxies should be more than 0. Hence, the mean of the random error must be set to 0 to fulfill the first rule, and the standard deviation must be set to the square of the mean value of the sample grids divided by 6, because 99% of the values which obey the normal distribution will be located in the range from mean ±3 standard deviation. This kept the minimum value of the random error as 0, to fulfill the second rule. Subsequently, proxies with new random errors were obtained and gridded to produce the new CO₂ emissions map. The new gridded CO₂ emissions maps were compared to the original gridded maps, and the difference was estimated as follows:

$$\left\{\begin{array}{c}\begin{array}{c}Weigh{t}_{i,l}{}^{’}=\overline{Weigh{t}_{20\%,l}}+ {\xi }_i, {\xi }_i~N\left(0,\frac{\overline{Weigh{t}_{20\%,l}{}^2}}{6} \right)\\ differenc{e}_{i,l}=simulatio{n}_{i,l}- Gri{d}_{i,l}\end{array}\\ difference=\frac{{{\displaystyle \sum }}_{i=1}^{n}Differenc{e}_{i,\mathrm{l}}}{C_l}=\frac{{{\displaystyle \sum }}_{i=1}^{.n}\left|Gri{d}_{i,l}-simulatio{n}_{i,l}\right|}{C_l}\times 100\%\end{array}\right..$$

(4)

where Weight_i,l’ is the new value of spatial proxies set to a random error of the lth sector, $\overline{Weigh{t}_{20\%,l}}$ is the mean of the sampled grid of the spatial proxies from lth sectors, and ξ_i is the random error. difference_i,l is the absolute difference between the real and simulated emissions for the lth sector on the ith grid (l = 1, 2, 3, and 4, which represent residential, industrial, transport, and total emissions, respectively), C_l is the total CO₂ emissions from the lth sector, and Grid_i,l and simulation_i,l represent the total and simulated CO₂ emissions for the lth sector on the ith grid.

3. Results

3.1. High-Resolution Gridded CO₂ Emissions Map Patterns

3.1.1. High-Resolution Gridded CO₂ Emissions Map

In this study, we used the gridded CO₂ emissions maps constructed in previous studies [26,27,28]. Please see Tables S1–S3 (Supplementary Materials) for the composition and the heterogeneity analysis of the maps. The annual CO₂ emissions per capita of Jinjiang city was 8.12 t in 2013, and 10,000 yuan of Jinjiang’s GDP caused 1.219 t of CO₂ emissions. The total emission patterns of the 30 m and 500 m resolution CO₂ emissions maps were similar (Figure 2a,e). The highest CO₂ emission grids were located in the urban residential areas in the northern region, and in the industrial areas in the central region. Rural areas had low values of CO₂ emissions. Residential emissions were concentrated in the northern built-up area, with high values of CO₂ emissions being observed on the 500 m resolution map in this built-up area (Figure 2b,f). The distribution of industrial emissions was wider and more fragmented on the 500 m resolution map than on the 30 m resolution map (Figure 2c,g). Transport CO₂ emissions were clearly distributed along roads on both maps (Figure 2d,h).

3.1.2. Spatial Variation of CO₂ Emissions

We constructed semi-variograms of CO₂ emissions to understand their spatial variation [32] (Figure 3); detailed information of the method is shown in the Supplementary Materials. The well-known geo-statistical concepts of nugget and sill represent the values at which the semi-variogram almost intercepts the y-value, and the values at which the model first flattens out, respectively [33]. The greater the sum of nugget and sill, the larger the spatial variation of CO₂ emissions [34]. The results indicate that the spatial variation of emissions was much larger at 30 m resolution than at 500 m resolution. At 30 m and 500 m resolution, the largest spatial variations were both recorded for transport emissions among the different sectors (Figure 3d,h).

3.2. Propagated Uncertainties Caused by Gridded Model

3.2.1. Results of the Uncertainty Estimation, Using U1

The logarithms of the mean and standard deviation of spatial proxies for CO₂ emissions were obtained by fitting the PDFs (Figure 4). Smaller KS test values indicated better PDF fits. The KS test values of different maps and sectors were all less than 0.3, and were considered as good fits. The fitting performance of the 500 m resolution map was better than that of the 30 m resolution map. Based on the PDF fitting, we drew the curves of simulated and real CDFs for different resolutions (Figure 4). The CDFs represents the number (i.e., cumulative probability density) of each grid value occupied with the total number of grids. The CDF curve is an ascending order by x value. Hence, the cumulative probability density of a larger x value is almost 1. The uncertainty was the area between the two curves (gray area in Figure 4). The uncertainty of any pixel value was the difference between the y coordinates of the two curves when the pixel value was the corresponding x value. Comparison of the CDF curves of simulated emissions with actual emissions (Figure 4) showed that the uncertainty of the residential emissions was the smallest, while that of transport emissions was slightly larger, and that of industrial emissions was the largest. Overall, the uncertainty in the 500 m resolution map was smaller, and tended to decrease as CO₂ emissions increased. In comparison, the uncertainty in 30 m CO₂ emissions showed no significant change as CO₂ emissions increased. When comparing the different sectors, there was less uncertainty regarding residential emissions in both the 30 m and 500 m resolution maps, which was represented by the coincident two color curves (Figure 4b,f). Uncertainty in industrial emissions was very large (up to nearly 80%) at 20,000 t/pixel, but it gradually decreased, and stayed within 10% beyond 20,000 t. The uncertainty on the 30 m map was almost always maintained at about 40% (Figure 4c,g). When the transport discharge pixel value was between 300 and 600 t, the uncertainty was small, but then gradually decreased with increasing pixel values on the 500 m map. There was always about 30% uncertainty in the 30 m map (Figure 4d,h).

The uncertainty map for the gridded model was generated by a look-up table (Figure 5). Among the uncertainty maps of total and industrial CO₂ emissions (Figure 5a,c,e,g), the number of grids whose uncertainty reached 100% accounted for more than 80% of the area. However, most of the corresponding CO₂ emissions of these grids were equal to zero. Most grids below 100% located in the residential and industrial areas of Jinjiang city. For the uncertainty maps of residential CO₂ emissions (Figure 5b,f), almost half of the grids of the 30 m resolution uncertainty maps have low uncertainty (i.e., the average value of the uncertainty grids was 0.969%). With the other resolution map, most uncertainty grids were less than 25%. For the uncertainty maps of transport CO₂ emissions (Figure 5d,h), most non-100% uncertainty grids were distributed along the roads. Compared with the 500 m resolution uncertainty map, the 30 m resolution uncertainty map had more 100% uncertainty grids. Furthermore, other grids of the 30 m and 500 m resolution uncertainty maps were below 50%.

3.2.2. Results of the Uncertainty Estimation Using U2

The statistical terms and population distributions (sums of the grid values) were obtained based on the bootstrap sampling method described in the methods section. Then, the uncertainties and CIs were calculated (Figure 6,

Table 1. Uncertainty in CO₂ emissions (x 10³ t) from different sectors from the 30 m and 500 m resolution maps.

Source	95 CI (Range %)		Mean		Uncertainty
	30 m	500 m	30 m	500 m	30 m	500 m
Residential	2304 to 2970 (−11% to +14%)	2401 to 2799 (−8% to +8%)	2596	2596	12.84%	7.68%
Industrial	9744 to 18,170 (−29% to +32%)	11,917 to 15,591 (−13% to +14%)	13,759	13,731	30.62%	13.38%
Transport	240 to 357 (−19% to +20%)	274 to 323 (−8% to +8%)	297	298	19.67%	8.19%
Total	12,289 to 21,498 (−26% to +29%)	14,593 to 18,714 (−12% to +13%)	16,653	16,624	27.65%	12.40%

Notes: the uncertainty is equal to the mean value of CO₂ emissions for the selected grid in the bootstrap sampling method divided by the 95% CIs width of CO₂ emissions at the selected grid in the bootstrap sampling method. CI, confidence interval.

), the smallest uncertainty was obtained for residential emissions (12.84% and 7.68%), while CO₂ emissions from industry remained the main source of uncertainty. Transport emissions were relatively small, with a small uncertainty. The uncertainty of the 500 m resolution map was far less than that of the 30 m resolution map. Since the range of 95% CIs on the 500 m resolution map was much smaller than that on the 30 m resolution map, the estimated overall CO₂ emissions of the 500 m resolution map were also closer to the calculated value. The relationship between the overall uncertainty of emissions and the uncertainty of the three sectors was not a simple linear superposition, with it presenting a nonlinear change, because the overall uncertainty was less than the sum of the uncertainty of the three sectors. The estimated total emissions were similar than the normal distribution based on 10,000 simulations (Figure 6). Only the estimates of residential emissions in the 30 m resolution map showed a right-skewed distribution (Figure 6b).

3.3. Propagated Uncertainties Caused by Input

3.3.1. Results of the Uncertainty Estimation Using U3

The constructed uncertainty maps for activity levels showed that the mean value of uncertainties for the grids of the 30 m resolution map was 37.30%, while that of the grids of the 500 m resolution map was 33.10%. The 95% CIs and uncertainty for residential, industrial, transport, and total activity levels were 1573 to 3594 (−39.13% ~ 39.08%), 11,005 to 16,424 (−19.81% ~ 19.67%), 205 to 392 (−31.25% ~ 31.59%), and 13,697 to 19,490 (−17.52% ~ 17.36%) thousand tons, and 39.11%, 19.74%, 31.42%, and 17.44%, respectively. The 95% CIs of the activity levels of the different emission types were within ±2σ. The 95% CIs of the total emissions were within ±20%, and had the lowest uncertainty compared to all of the sectors. The number of pixels with high uncertainty was greater in the 30 m resolution maps compared to the 500 m resolution maps (Figure 7). The uncertainties in the middle industrial areas and the northeastern built-up areas were smaller in both the 30 m and 500 m resolution maps (uncertainty was below 30%). Large areas with low uncertainty in the northern part of the 500 m resolution map had high uncertainty in the 30 m resolution map. Low uncertainty grids were scattered and fragmented in the 30 m resolution map (Figure 7a), but presented a continuous distribution in the 500 m resolution map (Figure 7b).

3.3.2. Results of the Uncertainty Estimation Using U4

U4 represents the uncertainty maps made with spatial proxies. The value of grids in the 30 m resolution map, which were regarded as the wrong grids, were generated by the product of average value of selected grids and a random error for residential, industrial, and transport emissions, were 19.85%, 20.00%, and 19.86%, respectively, while the total error was 33.67%. The errors in the spatial proxies of the distribution of residential, industrial, and transport emissions were 19.98%, 19.98%, and 19.98%, respectively, in the 500 m resolution map, while the total error was 24.35%. The overall error was greater than any initial error, and was less than the sum of the errors for the three sectors. The overall error presented a nonlinear change. The 500 m resolution map had less uncertainty than the 30 m resolution map. When there was a random error in spatial proxies, most areas overestimated emissions, with the errors of the urban areas and industrial areas with high emission intensity being large (Figure 8).

4. Discussion

4.1. Uncertainties

The process of mapping CO₂ emissions is characterized by various uncertainties. Gurney [35] identified challenges in quantifying uncertainties in CO₂ emissions, and suggested solutions, such as sensitivity testing of key inputs and parameters that are fully propagated through the estimation system. In this study, an analytic workflow inspired by the concept suggested by Gurney [35] was proposed to determine uncertainties at the urban scale.

Differences between estimations and actual emissions were generally used to measure the overall uncertainty in the estimation of emissions. Ou et al. [36] used relative error, mean relative error, and root mean square error to evaluate differences between actual and estimated data. This study used the CDF, which provided a description of the grid distribution features and allowed a comparison of the differences. The results indicated less uncertainty for residential emissions, which was in contrast to the uncertainties obtained from global mapping products. Andres et al. [21] analyzed uncertainty in three inputs of the CDIAC data product, and attributed 51% of the uncertainty to population data. Therefore, more accurate and detailed data were used to construct the population density gridded maps and reduce the uncertainties in this study. The results suggested that industrial emissions were the main uncertainty source in estimating urban CO₂ emissions.

Several studies have investigated uncertainties caused by total CO₂ emission estimations that are known to generate fundamental uncertainties [18,37]. Inventory accounting generates the quantities of total emissions that are eventually mapped. Since the uncertainties in activity levels and emission factors vary by sector and country, respectively, this study set the a priori uncertainty with normally distributed activity levels referenced by the paper of Zhao et al. (2012) [19]. Moreover, this study reported a greater, 95% CIs for total emissions than [19] did (−9 to 11%). The uncertainty estimated in this study (17.4%) was less than that determined from the CDIAC data for China (17.5%) [21]. These results could be attributed to the following reasons. First, this study focused on the propagation of uncertainty in activity levels during the gridding process. In contrast, Zhao et al. (2012) quantified uncertainty in the total CO₂ emissions estimation of China through providing a detailed inventory while without quantifying uncertainties for the gridded process. Second, the trends for uncertainties in the total emissions and in the sums of emissions from different sectors were similar to those reported by Zhao et al. (2012) [19]. The uncertainty caused by the total emission estimations was significantly reduced by calculating the detailed inventories from different sectors, likely due to the “compensation of error” phenomenon that was realized through Monte Carlo simulations [19]. Finally, in contrast to Zhao et al. (2012), the higher uncertainty detected in this study could be attributed to the lack of detailed source categories. However, a comparison of results from the current study and the CDIAC indicated acceptable uncertainty levels.

Geospatial datasets were primarily used to obtain more accurate spatial proxies and reduce uncertainties in the gridded process. Uncertainties caused by spatial proxies also varied across different sectors. For instance, uncertainties in residential emissions were mainly attributed to errors in the generation of high-resolution population grids [21]. However, high population fitting accuracy at the street and town scales was generated in this study, and such spatial proxies of population were used to map residential emissions. Additionally, uncertainties in industrial emissions were mainly driven by the accuracy of industrial point-source data identification and the extent to which night time light intensity represented emissions intensity (saturation of bright lights [4]). Raupach et al. [38] indicated that correcting the saturation error of nighttime light imagery increased data values by 1.15–1.23 times, leading to greater uncertainty. Uncertainties in transport emissions existed, since road areas do not fully represent the intensity of CO₂ emissions from traffic. Previous studies used traffic congestion conditions to calculate this parameter [39,40,41]. However, such data are difficult to acquire. Therefore, road area is generally used as a proxy in this research. Although road area is correlated with traffic jams and vehicle speed, it might not fully reflect traffic conditions [42].

In this study, data from the three sectors were ultimately spatially superimposed during the gridded process and represented the process of assigning spatial proxies. The U4 set 20% of the random errors in spatial proxies. Therefore, uncertainties caused by the total emission estimations were driven by the superimposition of interference from different sectors, were significantly less than 60%, and indicated a nonlinear increase. These results are supported by previous research on error transfer in cellular automaton models [43]. Yeh et al. [43] proposed that the errors in data sources did not propagate entirely to the final results. This means that the errors may be reduced after model processing, because their effects significantly decline during simulations. This argument is based on the assumption that the neighborhood calculation principle smooths errors. Therefore, the random error of total emissions was reduced by the compensation-of-error phenomenon, which is attributed to overlaying emissions from different sectors. Although the total quantity of emissions was stable for the whole city, errors increased at finer spatial resolutions. This finding is consistent with the results of Dong et al. [44], who mapped uncertainties in land use impacts on ecosystem services. Therefore, the impacts of spatial proxy errors on uncertainty were lower in coarse-resolution maps. This might be because fine-resolution maps have larger spatial variation in CO₂ emissions than coarser-resolution maps. Dale et al. [45] attributed uncertainty in their model to an inadequate representation of spatial variability in natural elements. Spatial scale has been shown to significantly impact the spatial variation and distribution patterns of natural elements [46]. In this study, the fine-resolution map indicated higher spatial variation and larger uncertainty (Figure 3) than the coarse-resolution map. These results are supported by Ogle et al. [47], who reported that uncertainty was negatively correlated with spatial scale in their study on soil organic carbon. Cias et al. [48] analyzed uncertainty in CO₂ emissions at the sub-national scale for 25 member states of the European Union, and concluded that uncertainty significantly increased when the resolution was smaller than 200 km. Therefore, multi-scale research is required to establish critical thresholds for scale selection [49,50,51].

The analytic workflow for uncertainty propagation proposed in this study may aid decision-makers in designing suitable low-risk policies. Regions with high emissions values and low uncertainties (industrial sources with emissions greater than 20,000 t/pixel and residential areas) could be the primary targets for emission reductions (Figure 4c). Therefore, the northern urban living areas and central regions could reduce emissions by promoting the development of low-carbon communities, and implementing energy-conservation and emission-reduction technologies. Furthermore, this study contributed to the research field of uncertainty analysis of urban scale CO₂ emissions by proposing an analytic and reproducible workflow consisting of four submodules to investigate the uncertainty propagation during the gridding process caused by the gridded model and input. Previous studies used Monte Carlo simulation to estimate the uncertainty in CO₂ emissions caused by inputs [10,11,12,15,52,53], the remaining studies used inventory comparison analysis to evaluate this uncertainty [5,25]. However, few studies considered the uncertainties caused by both the gridded model and the input. Furthermore, inventory comparison analysis is impossible at the urban scale due to the lack of a detailed emission inventory. The studies developed different uncertainty metrics and tools such as PDFs, the bootstrap sampling method, fuzzy mathematics, the 95% Cis, and relative uncertainty [54,55,56,57], which could provide the reference for our workflow. Moreover, some studies which retrieved CO₂ emissions using atmospheric models used data assimilation technologies in an attempt to estimate and reduce the uncertainty in CO₂ emissions [53,58]. These are novel technologies that can estimate the uncertainty but require complex atmospheric models and high-quality flux monitoring data. Compared with previous related studies, involving the construction and uncertainty propagation analysis of gridded maps of CO₂ emissions, the workflow combined various methodologies of uncertainty analysis, including Monte Carlo simulation and the bootstrap sampling method, and enabled finding the completed uncertainties caused by the gridded model and input without detailed emission inventories, field measurements, or integration of complex atmospheric models.

4.2. Future Research Directions

Gridded CO₂ emissions maps are likely to play an important role in climate change mitigation and sustainable development. Therefore, three future directions in this regard are proposed in the following sections based on a review of previous studies.

Advances in earth observation satellites and related technologies have enabled accurate, high-resolution quantification of CO₂ emissions. Gridded CO₂ emissions maps can be developed based on data obtained from night time light and greenhouse gas monitoring satellites [59,60,61]. However, relationships between night time light intensity and CO₂ emissions are largely empirical [4,62]. Greenhouse gas monitoring satellites determine the actual concentrations of the CO₂ column at low resolutions [63,64,65]. These limitations can be addressed by advances in earth observation satellites and implementation of novel technologies, including the development of satellites such as the LuoJia 1-01 and TanSAT [66,67,68]. However, limitations in estimating uncertainties in satellite-retrieved products still exist. The U4 which was developed in this study regarded a priori error of spatial input data or the quality assessment of satellite products as the uncertainty to estimate the uncertainties of gridded CO₂ emissions maps. When we used remote sensing satellite products (such as MODIS-derived NDVI products) as input data, the satellite’s quality assessment files could be used to define a priori error, and to subsequently conduct uncertainty analyses [69].

The combination of carbon flux monitoring data and atmospheric models is essential for developing spatially explicit bottom-up CO₂ emissions models. The bottom-up approach can estimate CO₂ emissions with a high accuracy [22]. However, it requires high-quality and spatiotemporal resolution input data. For instance, the use of ¹⁴CO₂ measurements can significantly assist the estimation of emissions from biological and fossil carbon pools [70,71,72]. Carbon flux monitoring data include anthropogenic and biogenic CO₂ emissions [73]. However, long-term data on CO₂ emissions and carbon fluxes are limited. Therefore, long-term carbon flux monitoring networks need to be established [74]. Although the Hestia and Anthropogenic Carbon Emissions System products have been used to develop spatially explicit CO₂ emissions models, most data were obtained from the Indianapolis, Salt Lake City, Baltimore, and Los Angeles flux experiments [73]. However, such experiments are uncommon, and these studies coupled carbon flux monitoring data with atmospheric or chemistry transport models, such as the weather research and forecasting (WRF) model, the WRF model coupled with chemistry, the WRF Model coupled with the vegetation photosynthesis and respiration model, etc. [75,76]. Therefore, future studies (especially in China) should focus on different combinations of carbon flux monitoring data and atmospheric models.

Data assimilation technology allows the integration of satellite observation and model output data from atmospheric and oceanic models. Some studies have applied four-dimensional data assimilation of atmospheric CO₂ using atmospheric infrared sounder observations [77]. However, studies focusing on data assimilation for CO₂ emissions are limited because the coupled atmospheric model is difficult to construct and run. Moreover, there is a lack of satellite observations at a fine spatiotemporal resolution, and/or carbon flux monitoring data. The key components of data assimilation include observation data, uncertainty in observations, model outputs, uncertainty in model outputs, and data assimilation algorithms. Therefore, our workflow for uncertainty analysis is a powerful tool to aid data assimilation. Future research could focus on the integration of the present workflow with different data assimilation algorithms.

5. Conclusions

In this study, an analytic workflow was established to conduct uncertainty analyses. The workflow included four submodules which were used to evaluate the uncertainties which were propagated during the gridding process. The first and second submodules were designed to investigate the uncertainties of gridded CO₂ emissions maps caused by the gridded model. The third and fourth submodules were used to analyze the uncertainties of the gridded CO₂ emissions maps caused by the input. The results of our study demonstrated that different uncertainties were associated with the gridded process according to the statistical graph and gridded maps. Furthermore, we found that the uncertainties in the coarse-resolution map had a lower contribution to the uncertainties caused by the total emission estimations and spatial proxies. This may have been due to the fact that the fine-resolution map had a larger spatial variation in CO₂ emissions. Furthermore, a nonlinear mapping relationship between the sum of the uncertainties of the gridded CO₂ emissions maps among the different sectors and the final estimated uncertainties was found when the uncertainties of different sectors were spatially superimposed. This nonlinear relationship was attributed to the “compensation of error” phenomenon. That is, the overestimated and underestimated values among different sectors were offset at the same grid. In conclusion, this study provides the corresponding uncertainties of baseline data for developing low-carbon cities, as well as a general workflow for analyzing uncertainties in gridded CO₂ emissions maps.