Anthropogenic CO₂ emission estimates in the Tokyo metropolitan area from ground-based CO₂ column observations

Urban areas are responsible for more than 40 % of global energy-related carbon dioxide (CO₂) emissions. The Tokyo metropolitan area (TMA), Japan, one of the most populated regions in the world, includes various emission sources, such as thermal power plants, automobile traffic, and residential facilities. In order to infer a top–down emission estimate, we conducted an intensive field campaign in the TMA from February to April 2016 to measure column-averaged dry-air mole fractions of CO₂ (XCO₂) with three ground-based Fourier transform spectrometers (one IFS 125HR and two EM27/SUN spectrometers). At two urban sites (Saitama and Sodegaura), measured XCO₂ values were generally larger than those at a rural site (Tsukuba) by up to 9.5 ppm, and average diurnal variations increased toward evening. To simulate the XCO₂ enhancement (ΔXCO₂) resulting from emissions at each observation site, we used the Stochastic Time-Inverted Lagrangian Transport (STILT) model driven by meteorological fields at a horizontal resolution of ∼1 km from the Weather Research and Forecasting (WRF) model, which was coupled with anthropogenic (large point source and area source) CO₂ emissions and biogenic fluxes. Although some of the diurnal variation of ΔXCO₂ was not reproduced and plumes from nearby large point sources were not captured, primarily because of a transport modeling error, the WRF–STILT simulations using prior fluxes were generally in good agreement with the observations (mean bias, 0.30 ppm; standard deviation, 1.31 ppm). By combining observations with high-resolution modeling, we developed an urban-scale inversion system in which spatially resolved CO₂ emission fluxes at >3 km resolution and a scaling factor of large point source emissions were estimated on a monthly basis by using Bayesian inference. The XCO₂ simulation results from the posterior CO₂ fluxes were improved (mean bias, −0.03 ppm; standard deviation, 1.21 ppm). The prior and posterior total CO₂ emissions in the TMA are 1.026 ± 0.116 and 1.037 ± 0.054 Mt-CO₂ d⁻¹ at the 95 % confidence level, respectively. The posterior total CO₂ emissions agreed with emission inventories within the posterior uncertainty, demonstrating that the EM27/SUN spectrometer data can constrain urban-scale monthly CO₂ emissions.

1 Introduction

The steady increase of atmospheric greenhouse gas (GHG) concentrations has accelerated recent climate change. Although urban areas account for only 2 %–3 % of the global land surface, approximately 44 % of energy-related carbon dioxide (CO₂) emissions come directly from urban areas (Seto et al., 2014). Urban areas are thus a main target for emission reductions, and many cities around the world have committed to reducing their GHG emissions through both the C40 Cities Climate Leadership Group (https://www.c40.org/, last access: 25 July 2023) and city-specific programs. To support urban emission reduction strategies, a variety of efforts are currently underway to build emission inventories with high spatial and temporal resolution. These inventories estimate GHG emissions by using a bottom–up approach, in which the total emissions from each source category are calculated by multiplying activity data (e.g., fuel consumption, traffic counts, housing statistics) by emission factors indicating GHG emissions per unit of activity. Because such detailed inventories have been developed only for certain cities, another type of global emissions database has been developed that relies on spatial proxies (e.g., night lights, population) to downscale total emissions at national or sub-national scales. Gurney et al. (2019), however, have reported large discrepancies between downscaled and bottom–up fossil fuel CO₂ emissions at the urban scale, mainly due to large point sources and road traffic.

Independent verification of these emissions datasets is highly desirable, and atmospheric observations are increasingly being used for this purpose. Emissions can be estimated at fine scale from atmospheric observations of GHG concentrations by using high-resolution atmospheric transport models in a top–down inversion approach. Networks of in situ GHG observation stations providing both operational observations and emission estimates have been established in several megacities (e.g., McKain et al., 2012; Lian et al., 2022; Lauvaux et al., 2016; Sargent et al., 2018). In addition, emission estimates have been obtained by conducting aircraft-based measurement campaigns over megacities (Ahn et al., 2020; Pitt et al., 2022) and by using laser absorption spectroscopy for 2-D tomographic measurements (Lian et al., 2019). For the Tokyo metropolitan area (TMA), Japan, the most populous metropolitan area in the world, Pisso et al. (2019) estimated anthropogenic CO₂ fluxes during the winters of 2005–2009 at a spatial resolution of 20 km × 20 km from in situ measurements obtained by tower- or ground-based instruments (Tsukuba, Kisai, and Dodaira) and commercial aircraft-based instruments (over Narita) by an inverse analysis with a Lagrangian transport model. Recently, the number of tower- and ground-based observation sites in the TMA have been expanded, and additional atmospheric components and isotopes are being measured (Sugawara et al., 2021).

Compared with surface point measurements, total column measurements are less sensitive to changes in the height of the planetary boundary layer (PBL) (Olsen and Randerson, 2004; McKain et al., 2012). Therefore, column measurements help to both mitigate PBL height errors in an atmospheric inversion system (Gerbig et al., 2008) and disentangle the effects of atmospheric mixing from the exchange of carbon between the surface and the atmosphere (Wunch et al., 2011). In addition, column data obtained in urban areas include information about emissions over a broader spatial domain than surface point data. Babenhauserheide et al. (2020) have estimated CO₂ emissions from Tokyo by conducting a statistical analysis of long-term measurements of column-averaged dry-air CO₂ mole fractions (XCO₂) obtained with a ground-based high-resolution Fourier transform spectrometer (FTS, IFS 125HR, Bruker Optics) located at Tsukuba, about 50 km from Tokyo, together with wind data obtained from operational radiosonde observations. GHG emissions from several megacities have been characterized by field campaigns conducted with multiple portable FTSs (EM27/SUN, Bruker Optics) (e.g., Hase et al., 2015; Makarova et al., 2021). Comparison of observed XCO₂ and column-averaged dry-air methane mole fraction (XCH₄) values with simulation results obtained by high-resolution transport modeling demonstrated that the simulations could capture XCO₂ and XCH₄ gradients between upwind and downwind sites (Vogel et al., 2019; Zhao et al., 2019, 2023). Furthermore, emissions from megacities have been estimated from XCO₂ and XCH₄ data with Lagrangian transport models (Ionov et al., 2021; Jones et al., 2021; Hedelius et al., 2018). Cusworth et al. (2020) have estimated spatially resolved CH₄ emissions in the Los Angeles Basin at a resolution of 3 km × 3 km from operational surface and column data. Meanwhile, permanent city observatories with the EM27/SUN spectrometers are emerging (e.g., Dietrich et al., 2021).

In the present study, we performed an observation campaign using two EM27/SUN FTSs and the Tsukuba IFS 125HR FTS to constrain CO₂ emissions around Tokyo during late winter and early spring, when the proportion of clear days is high and the contribution of the biogenic flux to CO₂ fluctuations is minor. We constructed CO₂ emission inventories with more accurate information on both the locations and emissions of large point sources. Anthropogenic CO₂ emissions from area sources and large point sources were estimated separately using this inventory as the prior. In addition, the area source emission estimates with higher spatial resolution enable verification of the emissions reported by each administrative division. To simulate atmospheric transport at high spatiotemporal resolution, we used a Lagrangian transport model, the Stochastic Time-Inverted Linear Transport (STILT) model, driven by the Weather Research and Forecasting (WRF) model (Lin et al., 2003; Nehrkorn et al., 2010). We estimated spatially resolved anthropogenic CO₂ emissions by an inverse analysis and then evaluated the estimated monthly CO₂ emissions against inventory-based fossil fuel CO₂ emissions in the TMA. In Sect. 2, we describe the XCO₂ measurements by ground-based FTSs at three observation sites in the TMA. Section 3 presents a methodology for modeling XCO₂ enhancements at the observation sites using the atmospheric transport model and prior emission data, and a framework for estimating anthropogenic CO₂ emissions through Bayesian inference. In Sect. 4, we show the results of the XCO₂ measurements and simulations and discuss the posterior emission estimates.

Figure 1 Locations of the two EM27/SUN observation sites (Saitama University and Sodegaura City Hall) and Tsukuba TCCON site (yellow circles). Also shown are the AMeDAS stations (red circles) used for the comparison with wind data from the WRF simulation and the ERA5 reanalysis data. The calculation of the footprint by WRF–STILT was performed for the entire region displayed in this figure. The elevation data are from the Global Bathymetry and Topography at 15 arcsec (SRTM15 + V2.1) (Tozer et al., 2019). The upper-right figure shows the location of the study area relative to Japan as a whole.

2 Measurements during the observation campaign in the TMA

An observation campaign with portable FTSs and the Tsukuba FTS (the 2016 Tokyo campaign) was conducted in the TMA from February to April 2016. Here, the TMA is defined as a rectangular region that includes the urban areas of Kanagawa, Chiba, Saitama, Ibaraki, Tochigi, and Gunma prefectures as well as the Tokyo metropolis (Fig. 1). The United Nations reports that “Tokyo”, with approximately 38 million inhabitants, is the world's most populous area and accounts for 30.1 % of the total population of Japan (United Nations, 2019), although the boundaries used to define Tokyo by the United Nations are not the same as those used here. On the one hand, the city-center of the TMA, primarily because of its intense economic activity and high population density, is a strong source of anthropogenic CO₂ emissions (Saito et al., 2022), and there are a multitude of point sources with large emissions (large point sources), such as power plants and steel plants, along the shores of Tokyo Bay. On the other hand, the TMA is surrounded by evergreen broadleaf and deciduous broadleaf forests, which contribute to temporal variations in biogenic CO₂ fluxes. The 2016 Tokyo campaign was conducted from late winter to early spring, when biological activity within the TMA was dormant; thus, it was most likely responsible for much smaller changes in CO₂ concentrations than the anthropogenic activity.

We used measurement data from the ground-based high-resolution FTS operated as part of the Total Carbon Column Observing Network (TCCON, Wunch et al., 2011) at the National Institute for Environmental Studies (NIES) (35.0513^∘ N, 140.1215^∘ E, 31 m a.s.l. (above sea level)) in Tsukuba (Morino et al., 2018). Additionally, we used data from two EM27/SUN spectrometers (SN38 and SN44) throughout the campaign period and a third EM27/SUN spectrometer (SN63) beginning in the middle of the campaign. Considering the prevailing wind direction in this winter/early spring (i.e., northwesterly) and proximity to anthropogenic emission sources, we deployed the SN38 EM27/SUN spectrometer at Saitama University (35.8636^∘ N, 139.6081^∘ E, 28 m a.s.l.) and the SN44 EM27/SUN spectrometer at Sodegaura City Hall (35.4297^∘ N, 139.9545^∘ E, 14 m a.s.l.) from 16 February 2016 to 6 April 2016 (Fig. 1). The SN63 EM27/SUN spectrometer began operation at Tsukuba on 3 March 2016 and has since been continuously operated as part of the Collaborative Carbon Column Observing Network (COCCON, Frey et al., 2019). Before and after the observations at the three sites, we conducted side-by-side measurements with the EM27/SUN instruments and the TCCON instrument at NIES for approximately 1 week each time (3–10 February and 11–20 April).

The EM27/SUN instrument measures direct solar absorption spectra from 4000 to 11 000 cm⁻¹ in the short-wavelength infrared (SWIR) region (Gisi et al., 2012). At the time of the 2016 Tokyo campaign, the participating EM27/SUN spectrometers were equipped with only one InGaAs detector, and column amounts of CO₂, CH₄, water vapor, and oxygen were obtained from the SWIR spectra. The spectral resolution was 0.5 cm⁻¹, which corresponds to an optical path difference of 1.8 cm. Interferograms were continuously acquired every 6 s, and every 10 interferograms (five each from the forward and backward scans) were averaged and recorded (i.e., a sampling interval of about 1 min). If the weather permitted, EM27/SUN measurements at Saitama were made from sunrise to sunset, whereas the measurements at Sodegaura were made between approximately 09:00 Japan standard time (JST) and sunset, reflecting the office hours of the city hall. For retrieval processing, we used GGG2014 software, which is also used for processing the TCCON spectra (Wunch et al., 2015).

The EM27/SUN data were averaged in 15 min bins. Chen et al. (2016) derive an optimal integration time of 10–20 min, based on the Allan variance of two sets of EM27/SUN data from side-by-side measurements. However, they used a shorter integration time of 5 min to derive the EM27/SUN differences between upwind and downwind of local emission sources. In the present study, we found that it is difficult for the XCO₂ simulation to accurately reproduce the times at which point source plumes are observed (Sect. 4.2), and a comparison of the simulations and observations at short time intervals is not beneficial. Thus, we adopted an integration time of 15 min for the EM27/SUN data.

To ensure consistency among the instruments, we determined correction factors for XCO₂ and XCH₄ values based on the side-by-side measurements performed at NIES before and after the field campaign. The SN63 EM27/SUN data, which were bias-corrected against coincident aircraft measurements (Ohyama et al., 2020), were used as reference data. The instrumental line shape of the SN44 EM27/SUN deviated greatly from that of an ideal instrument because of its imperfect alignment (Frey et al., 2019); therefore, airmass-dependent correction factors (ADCFs) were derived in addition to airmass-independent correction factors (AICFs). The procedure for determining these correction factors is described in detail by Ohyama et al. (2021). The resulting correction factors C₀ (AICF) and C₁ (ADCF) in Eq. (1) of Ohyama et al. (2021) were 1.0028 and 0.0096, respectively, for the SN44 EM27/SUN XCO₂ data. The C₀ value of the SN38 EM27/SUN XCO₂ data was 1.0006. The C₀ and C₁ values of SN44 EM27/SUN XCH₄ data were 1.0101 and 0.0021, respectively, and the C₀ value of the SN38 EM27/SUN XCH₄ data was 1.0017. Comparisons of the bias-corrected EM27/SUN data with the SN63 EM27/SUN data are shown in Fig. S1 in the Supplement. Each of the EM27/SUN data points was averaged per 15 min bin. After the bias correction, the SN38 and SN44 EM27/SUN XCO₂ data differed from the SN63 EM27/SUN XCO₂ data by (mean ± 1σ) −0.01 ± 0.17 and 0.06 ± 0.16 ppm, respectively, and the SN38 and SN44 EM27/SUN XCH₄ data differed from the SN63 EM27/SUN XCH₄ data by −0.09 ± 0.97 and 0.39 ± 0.88 ppb, respectively. In the present study, the Tsukuba TCCON data were also scaled to be consistent with the SN63 EM27/SUN data by using C₀ values of 0.9977 for XCO₂ and 0.9985 for XCH₄. Figure 2 shows time series of the bias-corrected XCO₂ and XCH₄ data observed with the four spectrometers during the campaign period, including the side-by-side measurements at the Tsukuba site.

Figure 2 Time series of XCO₂ and XCH₄ during the observation campaign, including side-by-side measurements conducted at Tsukuba. The dashed vertical lines show the dates when the field observations began and ended at the three sites.

3 Methodology for the CO₂ emission estimation

3.1 Lagrangian transport model

We used the STILT model (Lin et al., 2003; Fasoli et al., 2018) coupled with meteorological fields from the WRF model (hereafter WRF–STILT; i.e., STILT driven by meteorological fields from WRF) to simulate atmospheric transport as required for inversion of the CO₂ emission data. STILT calculates the trajectory of particles from a “receptor” location and generates a footprint that represents the sensitivity of the CO₂ mole fraction to be measured at the receptor location to upwind emissions. This footprint (concentration per unit flux; ppm (mol m⁻² s⁻¹)⁻¹)) corresponds to the Jacobian matrix used for inverse analysis of CO₂ emissions. We used WRF meteorological fields generated at a horizontal resolution of approximately 1 km to drive STILT (described in detail in Sect. 3.2). We ran the model at 14 receptor levels (25, 50, 75, 100, 150, 200, 300, 400, 600, 800, 1000, 1500, 2000, and 2800 m a.g.l.) over each observation site (Jones et al., 2021), and an ensemble of 1000 particles for each altitude was traced backwards in time for 24 h. We varied the location (latitude and longitude) of the receptor along the line-of-sight of the EM27/SUN pointing toward the sun to accord with the receptor level. The hourly footprints for each grid were computed by considering the PBL height and the residence time of particles traveling within the lower PBL (Lin et al., 2003). The footprint calculations were performed every 15 min from 09:00 to 17:00 JST at a spatial resolution of 30 arcsec (approximately 1 km) within 34.975–36.625^∘ N, 138.900–140.875^∘ E (Fig. 1).

The hourly footprints calculated over the STILT run time (24 h) at a given time were weighted by temporal correction factors of CO₂ emissions (described in Sect. 3.3) and aggregated in each grid cell. From the summed footprints at each altitude, we then calculated the pressure-weighted column-average footprint, taking account of the column-averaging kernel of the EM27/SUN spectrometer (Rodgers and Connor, 2003; Jones et al., 2021). The footprints generated by STILT were then multiplied with spatially resolved emission inventories for anthropogenic and biogenic fluxes separately to determine the spatially resolved contributions (in ppm) of the surrounding emission sources. The change in XCO₂ (ΔXCO₂) at each observation site was obtained by summing the contributions over all grid cells and serves for the forward modeling. We separated the anthropogenic flux into large point source and area source emissions, as described in Sect. 3.3. We thus considered three types of fluxes: large point source emissions, area source emissions, and the biogenic flux.

3.2 Meteorological fields from WRF model

To drive the STILT model, we used meteorological fields simulated using the Advanced Research WRF model (WRF–ARW version 3.9.1.1; Skamarock et al., 2008). As meteorological initial and lateral boundary conditions for the WRF simulation, we used grid point value (GPV) data produced by the mesoscale forecast model (MSM) of the Japan Meteorological Agency (JMA) (MSM–GPV; JMA, 2019). The MSM–GPV data have 17 vertical layers, including the surface layer with a spatial resolution of 0.0625^∘ × 0.0500^∘ and 16 pressure levels (from 1000 to 100 hPa) with a spatial resolution of 0.125^∘ × 0.100^∘. Although the MSM–GPV data provide 3-hourly forecast values, only the initial values of each forecasting cycle were used in this study. The initial and lateral boundary conditions of soil temperature and moisture were obtained from final operational global analysis and forecast data (GDAS/FNL) of the National Centers for Environmental Prediction (NCEP) with a spatial resolution of 0.25^∘ × 0.25^∘ (NCEP, 2015). Because the spatial resolution of the MSM–GPV data at the 16 pressure levels was ∼10 km, the WRF model was configured with two modeling domains (d01 and d02 in Fig. S2) with horizontal resolutions of 3 km (d01) and 1 km (d02), with one-way grid nesting. Domain 01 included not only the Kanto Plain but also the surrounding mountainous and marine areas, and domain 02 fully covered the TMA and was slightly larger than the domain used for the STILT simulations. We set 51 vertical levels extending from the surface up to 100 hPa. Land use information was taken from a dataset (veg_jstream) constructed by Japan's Study for Reference Air Quality Modeling (Chatani et al., 2018). To reduce computational effort, the WRF simulations were not run for the entire campaign period but for separate intervals of 2.5–5.5 consecutive days, which were determined so that they covered the EM27/SUN measurement days. Each simulation segment started at 12:00 UTC, and the first 12 h was considered spin-up time. Grid nudging toward the MSM–GPV data was applied to the wind field (u,v), temperature (t), and the water vapor mixing ratio (q) at all levels in domain 01 with a nudging coefficient of 3.0 × 10⁻⁴ s⁻¹ for each element. In domain 02, grid nudging of the wind field was applied at all levels with a nudging coefficient of 1.0 × 10⁻⁴ s⁻¹, whereas nudging was applied to the temperature and water vapor mixing ratio only at the levels above the simulated PBL with a nudging coefficient of 3.0 × 10⁻⁵ s⁻¹. The simulations for domains 01 and 02 were carried out with integration time steps of 15 and 5 s, with model outputs saved every 1 h and every 15 min, respectively. For use with the STILT model, the wind data for domain 02 were time-averaged over the output interval of the WRF model (Nehrkorn et al., 2010). Table 1 summarizes the model settings and physics options used for the reference inverse analysis.

Table 1 Physics and model options of the WRF data used in the reference inversion.

In previous studies using the WRF model, the physics options of the model were set according to the studied region and period as well as the WRF version. In this study, we sought to select optimal physics options especially for the PBL scheme, the cumulus parameterization scheme, and the land surface model, all of which significantly impact the WRF calculation (Díaz-Isaac et al., 2018), by comparing WRF (and STILT) simulation results obtained with different physics options against measurement data (see Sect. 4.2 for the STILT simulation). The Kain–Fritsch cumulus parameterization scheme (Kain, 2004) was applied only in domain 01 (Table 1). Because cumulus parameterization is valid only for coarse grid resolutions (typically >10 km), we also ran simulations without any cumulus parameterization scheme and found little difference in the simulation results obtained with and without a cumulus parameterization scheme. For the land surface model, we adopted the Rapid Update Cycle (RUC) model (Smirnova et al., 2016) because XCO₂ simulations using the RUC model reproduced the observations better than the other land surface models examined by Díaz-Isaac et al. (2018). We evaluated in detail the effect of different PBL schemes on the WRF simulation results because which PBL scheme is optimal depends on the location and season (e.g., Jeong et al., 2013). We compared the modeled wind fields with observational data from Automated Meteorological Data Acquisition System (AMeDAS) stations operated by the JMA (https://www.jma.go.jp/jma/en/Activities/amedas/amedas.html, last access: 25 July 2023). Because wind fields directly influence atmospheric transport patterns, it is of particular importance to assess the model performance of the wind fields. Here, we compared wind speed and wind direction in WRF simulations among three different PBL schemes, the Mellor–Yamada–Janjić (MYJ) scheme (Janjić, 1994), the Mellor–Yamada Nakanishi Niino Level 2.5 (MYNN25) scheme (Nakanishi and Niino et al., 2009), and the Yonsei University scheme (Hong et al., 2006) with topographic correction for surface winds (Jimenez and Dudhia, 2012) (YSU + topo), and we also compared the fifth-generation atmospheric reanalysis (ERA5) data at 0.25^∘ spatial resolution (Hersbach et al., 2023) with data from the five AMeDAS stations within the TMA (Fig. 1): Saitama (35.8761^∘ N, 139.5868^∘ E, 18 m a.s.l.), Tokyo (35.6916^∘ N, 139.7532^∘ E, 56 m a.s.l.), Haneda (35.5636^∘ N, 139.7896^∘ E, 16 m a.s.l.), Chiba (35.6028^∘ N, 140.1040^∘ E, 51 m a.s.l.), and Kisarazu (35.3623^∘ N, 139.9402^∘ E, 68 m a.s.l.). Comparison of wind speed and wind direction between the models and observations (see Fig. S3 for the Tokyo site) showed that the model data could generally reproduce the observed temporal variations in the wind fields and demonstrated the capability of the model to simulate reasonable meteorological fields for driving the transport of trace gases. Tables 2 and 3 summarize the mean differences (biases) between the models and observations and their standard deviations for wind speed and wind direction, respectively. The WRF MYJ scheme had the lowest bias in wind speed and the smallest standard deviation in wind direction, whereas the ERA5 results were the best for the standard deviation in wind speed and the bias in wind direction. Among the tested PBL schemes of the WRF model, the wind fields obtained with the MYJ scheme were optimal. XCO₂ simulations using these meteorological fields are assessed in Sect. 4.2.

Table 2 Mean differences and their standard deviations (1σ) in wind speed (m s⁻¹) between model data (three WRF simulations and ERA5 reanalysis data) and the observational data at five AMeDAS sites (model minus AMeDAS). Boldface indicates the best-case results among the models.

Table 3 Mean differences and their standard deviations (1σ) in wind direction (degrees) between model data (three WRF simulations and ERA5 reanalysis data) and the observational data at five AMeDAS sites (model minus AMeDAS). Boldface indicates the best-case results among the models.

Figure 3 Anthropogenic CO₂ emissions from the TMA in March 2016 in (a) the original ODIAC 2020b data and (b) the same data except that the locations and emission magnitudes of large point sources, such as power plants and manufacturing plants, were corrected based on the national emission inventory. Open diamonds denote the locations of large point sources, and the crosses in (b) denote the observation sites.

3.3 Anthropogenic and biogenic CO₂ fluxes

For the prior estimate of anthropogenic CO₂ emissions, we used the 2020b version of the Open-source Data Inventory for Anthropogenic CO₂ (ODIAC 2020b), which is a global high-resolution (30 arcsec) monthly fossil fuel CO₂ emissions database (Oda and Maksyutov, 2015; Oda et al., 2018). The high-resolution ODIAC emission map was created by spatially disaggregating the national CO₂ emissions using the large point source database and proxy data associated with emissions. The locations and magnitudes of large point source emissions in ODIAC are taken from the Carbon Monitoring for Action (CARMA) database (https://www.cgdev.org/topics/carbon-monitoring-action, last access: 25 July 2023). The rest of the national emissions (referred to as “area source emissions” because line source emissions are not explicitly included in ODIAC) are distributed based on the spatial patterns of night lights data collected by satellites. We pinpointed large point sources such as power plants and steel plants on the CO₂ emission map of ODIAC 2020b in March 2016 with the aid of high-resolution satellite imagery (Fig. 3a). The locations of large point sources in the ODIAC are not exact, most likely because of the large uncertainty of the CO₂ emission source information in the CARMA database (Gurney et al., 2019). We therefore customized the ODIAC data by a two-step process. First, grid cells with CO₂ emissions larger than a given threshold (>10⁴ ton of carbon per month) were replaced with the averaged value of the neighboring eight grid cells. We regarded these secondary emissions as area source emissions. Second, annual emissions from large point sources, which are available upon request from the Ministry of the Environment of Japan under the GHG Emissions Accounting, Reporting, and Disclosure System (https://ghg-santeikohyo.env.go.jp/, last access: 25 July 2023), were converted to monthly values and added to the area source emissions. The emission map corrected for the large point sources (referred to as the “LPS-corrected ODIAC”; Fig. 3b) was used as the prior estimate. Large point source and area source emissions comprised 37 % and 63 %, respectively, of the LPS-corrected ODIAC data. In addition to this correction, weekly and diurnal scaling factors from the Temporal Improvements for Modeling Emissions by Scaling (TIMES) model developed by Nassar et al. (2013) were applied to the ODIAC data to temporally downscale the monthly ODIAC product.

A bottom–up fossil fuel emission inventory in Japan with a spatial resolution of approximately 1 km × 1 km, the Multiscale Overlap Scheme for Analyzing national Inventory of anthropogenic CO₂ (MOSAIC), was used in a complementary manner (described in Sect. 4.3) (Saito et al., 2022). This emission inventory provides monthly data obtained by using Japanese government statistics for all socioeconomic activities of Japanese society. The inventory comprises fossil fuel CO₂ emissions from eight sectors: electricity generation, waste incineration, civil aviation, waterborne navigation, road transportation, industrial and commercial sources, residential sources, and agricultural machine use. The locations and emission magnitudes of the large point sources in the electricity generation sector were corrected using the same method applied to the ODIAC data. Note that this study used MOSAIC emission data for 2015, which were the only MOSAIC data available when the study was carried out.

To account for the influence of biogenic CO₂ on the observed ΔXCO₂ values, we used biogenic CO₂ fluxes calculated by a terrestrial biosphere model. Specifically, hourly net ecosystem exchange (NEE) data from the Vegetation Integrative SImulator for Trace gases (VISIT) model, referred to as “VISITc”, were adopted as the biogenic CO₂ flux data. The NEE data were combined with green vegetation fraction (GVF) data to downscale them. The NEE data reflect the CO₂ flux between the terrestrial biosphere and the atmosphere and are obtained as the difference between ecosystem respiration (RE) and gross primary productivity (GPP) in the VISITc data (RE–GPP). Whereas the initial VISIT products provided monthly fluxes with a 0.5^∘ spatial resolution (Ito and Inatomi, 2012), the VISITc products provide hourly resolved fluxes composited with meteorological input data from the Climate Forecast System Reanalysis (CFSR) versions 1 and 2 (Saha et al., 2010) and ERA-Interim (Dee et al., 2011). The VISITc model operates on the same grid as the CFSR data (i.e., approximately 0.31^∘ × 0.31^∘). The original VISIT model simulates the terrestrial biogeochemical cycle with a monthly resolution considering minor carbon flows, such as the effect of land-use change and fire disturbance, that directly affect variability in GPP and RE (Ito, 2019). These effects, however, were not adopted in VISITc, and therefore this study scaled the GPP and RE data derived from VISITc to those of original VISIT in every month and every grid. Furthermore, to better characterize the spatial distribution of biogenic CO₂ fluxes, we spatially downscaled the hourly VISITc NEE data using GVF data from the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor onboard the Suomi National Polar-orbiting Partnership satellite (VIIRS Global Green Vegetation Fraction). The GVF data are produced with a spatial resolution of approximately 4 km on a daily basis from the past 7 d of VIIRS observations (Ding and Zhu, 2018). The GVF parameter, which represents how much downward solar insolation is intercepted by the canopy, is used as a scaling parameter for the biogenic flux. Following the method of Ye et al. (2020), the original GVF data (Fig. 4a) were first averaged into the VISITc grid. Then, the VISITc NEE data (Fig. 4b) and the re-gridded GVF data were interpolated bilinearly into the 1 km × 1 km grid. The ratio of the original GVF to the interpolated GVF was multiplied by the interpolated NEE data to produce the downscaled NEE data (Fig. 4c). Finally, the downscaling process was conducted in a manner that ensured all original sums of the NEE data from the TMA were preserved following the downscaling. We note that the effective spatial resolution of the downscaled biogenic fluxes is about 4 km, although they were generated on a 1 km × 1 km grid.

Figure 4 (a) NEE data at 03:00 UTC on 23 March 2016 from the VISITc model and (b) GVF data during 20–26 March 2016. (c) NEE data downscaled using the GVF data.

3.4 Forward model

XCO₂ measurements at a given location are quantitatively related to the presumed surface CO₂ fluxes via the forward model H,

$$\begin{array}{}\text{(1)}& \mathit{y}=\mathit{H}\left(\mathit{x},\mathit{b}\right)+\mathit{\epsilon },\end{array}$$

where y denotes the measurement vector (n × 1), x is the state vector to be retrieved (m × 1), b is the vector consisting of fixed physical quantities, and ε is the measurement error vector. In the present study, the state vector x includes spatially resolved fluxes for the area source emissions and a scaling factor for the large point source emissions. Because the geolocations of large point sources are precisely prescribed (at the grid cell scale of ∼1 km), the emissions from the large point sources are treated separately from the area source emissions, which have large uncertainty in their spatial distribution. To simplify the inversion, the biogenic flux was allocated to the fixed vector b. Note that the contribution of biogenic flux to the simulated ΔXCO₂ was small compared to that of anthropogenic flux, and the differences among ΔXCO₂ calculated from four different biogenic flux products are also small (Sect. 4.2). The forward model simulates XCO₂ values at the urban sites (Saitama or Sodegaura) as follows:

$$\begin{array}{}\text{(2)}& \mathit{H}\left(\mathit{x},\mathit{b}\right)=\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{2}\mathrm{STILT}}^{\mathrm{urban}}\left(\mathit{x}\mathbf{,}\mathit{b}\right)+{\mathrm{XCO}}_{\mathrm{2}}^{\text{BG}}\left(\mathit{x}\mathbf{,}\mathit{b}\right),\end{array}$$

where $\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{2}\mathrm{STILT}}^{\mathrm{urban}}$ is the XCO₂ enhancement at the urban sites simulated by the pressure-weighted footprint and the surface fluxes, and ${\mathrm{XCO}}_{\mathrm{2}}^{\mathrm{BG}}$ is the background value. We calculated the ΔXCO₂ values as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{2}\mathrm{STILT}}^{\mathrm{urban}}\left(\mathit{x}\mathbf{,}\mathit{b}\right)={\mathit{F}}_{\text{aggr}}^{\mathrm{urban}}{\mathit{x}}_{\mathrm{area}}+{\mathit{F}}_{\mathrm{fine}}^{\mathrm{urban}}{\mathit{b}}_{\mathrm{point}}{x}_{\mathrm{point}}\\ & +{\mathit{F}}_{\mathrm{fine}}^{\mathrm{urban}}{\mathit{b}}_{\mathrm{bio}},\end{array}\end{array}$$

where F_fine and F_aggr are the original and the spatially aggregated footprints, respectively. x_area and x_point are the emission flux vector for area sources and the (scalar) scaling factor for large point sources, respectively. b_point and b_bio are the emission flux vectors for large point sources and biogenic sources, respectively. The ΔXCO₂ values resulting from the large point source emissions and biogenic fluxes were calculated from the original footprints with a spatial resolution of approximately 1 km × 1 km (0.0083^∘ × 0.0083^∘ ). For area source emissions, however, we re-gridded the original footprints to a spatial resolution of 0.025^∘ × 0.025^∘ to degrade the spatial resolution for the inverse analysis. First, the area source emissions were summed for each 0.025^∘ × 0.025^∘ grid cell. Then, individual footprints for the 0.025^∘ × 0.025^∘ grid were derived by dividing the sum of the nine XCO₂ contributions for the 0.0083^∘ × 0.0083^∘ grid by the emissions for the 0.025^∘ × 0.025^∘ grid.

We assumed that the XCO₂ values at Tsukuba approximately represent background air, as there are lower CO₂ emissions around Tsukuba (Fig. 3) and the XCO₂ values observed at Tsukuba were systematically lower than those at the other urban sites, which can be seen from the XCO₂ values in Fig. 2a. However, the observations at Tsukuba are not strictly background; they are affected by emissions from the TMA (Babenhauserheide et al., 2020). We therefore obtained the background XCO₂ values by subtracting the simulated ΔXCO₂ values at the Tsukuba site ($\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{2}\mathrm{STILT}}^{\mathrm{Tsukuba}})$ from the Tsukuba TCCON XCO₂ values (${\mathrm{XCO}}_{\mathrm{2}\mathrm{TCCON}}^{\mathrm{Tsukuba}})$:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathrm{XCO}}_{\mathrm{2}}^{\mathrm{BG}}\left(\mathit{x}\mathbf{,}\mathit{b}\right)={\mathrm{XCO}}_{\mathrm{2}\mathrm{TCCON}}^{\mathrm{Tsukuba}}-\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{2}\mathrm{STILT}}^{\mathrm{Tsukuba}}\left(\mathit{x}\mathbf{,}\mathit{b}\right)\\ & ={\mathrm{XCO}}_{\mathrm{2}\mathrm{TCCON}}^{\mathrm{Tsukuba}}-\left({\mathit{F}}_{\mathrm{aggr}}^{\mathrm{Tsukuba}}{\mathit{x}}_{\mathrm{area}}+{\mathit{F}}_{\mathrm{fine}}^{\mathrm{Tsukuba}}\right.\\ & \left.{\mathit{b}}_{\mathrm{point}}{x}_{\mathrm{point}}+{\mathit{F}}_{\mathrm{fine}}^{\mathrm{Tsukuba}}{\mathit{b}}_{\mathrm{bio}}\right).\end{array}\end{array}$$

The background values were presumed to be common to the three sites, given the proximity of the sites. When there were no data available from the Tsukuba site, CarbonTracker CT2019B XCO₂ data (CT2019B.xCO2) (Jacobson et al., 2020) were used for background values independent of local time. The average CarbonTracker data of two grids centered over the ocean east of the TMA (35.0^∘ N, 142.5^∘ E and 37.0^∘ N, 142.5^∘ E) were used. We note that from February to April 2016, the mean difference between the CarbonTracker data (13:30 LT) and the Tsukuba TCCON data (average of 12:00–15:00 LT) (CarbonTracker minus TCCON) was −0.18 ppm with a standard deviation of 0.72 ppm. In the following simulations and inverse analyses, only the TCCON data were used as the measurement data at Tsukuba, since the SN63 EM27/SUN measurements started in the middle of the campaign (as described in Sect. 2).

The simulations were performed when EM27/SUN measurements at Saitama or Sodegaura were available for more than 2 h per day and when two or more of the 15 min averaged data showed XCO₂ differences of at least 1 ppm between the urban (i.e., Saitama or Sodegaura) and Tsukuba TCCON site. At Saitama and Sodegaura, these conditions were satisfied at 15 and 11 d, respectively, between 16 February and 23 March 2016. In addition, the simulations were limited to the time period from 09:00 to 17:00 JST, when synchronous measurements at the three sites were made and the airmass-dependent variation in XCO₂ was moderate. Because the time period when the forward simulations and the subsequent inverse analysis were performed was approximately 1 month during February and March 2016, for the prior of the anthropogenic emissions we used the average of the ODIAC data in February and March 2016.

3.5 Inversion methodology

The anthropogenic CO₂ emission fluxes can be estimated from the observed XCO₂ values, together with their associated uncertainties, through an optimization procedure. For the area source emissions, we decided to optimize the logarithm of the emission flux, with the prior errors following a Gaussian distribution, because the area source emissions from each grid cell differ by a couple of orders of magnitude, and the optimization of area source emissions at linear scale might lead to unphysical negative posterior emissions. On the other hand, the scaling factor for the large point sources was optimized at linear scale. These anthropogenic CO₂ emissions (i.e., x) were estimated based on a Bayesian framework by minimizing the cost function, which consists of two terms related to the measurements and prior emissions:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& J\left(\mathit{x}\right)={\left(\mathit{y}-\mathit{H}\left(\mathit{x}\right)\right)}^T{\mathbf{S}}_{\mathit{\epsilon }}^{-\mathrm{1}}(\mathit{y}-\mathit{H}\left(\mathit{x}\right))\\ & +{\left(\mathit{x}-{\mathit{x}}_{\text{a}}\right)}^T{\mathbf{S}}_{\text{a}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}_{\text{a}}),\end{array}\end{array}$$

where x_a is the prior vector of x; S_ε is the model–observation mismatch covariance (or the measurement error) matrix (n × n); and S_a is the priori error covariance matrix (m×m). Because the state vector x includes the logarithm of area source emissions, the inverse problem is nonlinear, and an iterative approach was used to estimate the CO₂ emissions (Rodgers, 2000):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathit{x}}_{l+\mathrm{1}}={\mathit{x}}_l+{\left[{\mathbf{K}}^T{\mathbf{S}}_{\mathit{\epsilon }}^{-\mathrm{1}}\mathbf{K}+\left(\mathrm{1}+\mathit{\gamma }\right){\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right]}^{-\mathrm{1}}\\ & \left[{\mathbf{K}}^T{\mathbf{S}}_{\mathit{\epsilon }}^{-\mathrm{1}}\left(\mathit{y}-\mathit{H}\left({\mathit{x}}_l\right)\right)-{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}({\mathit{x}}_l-{\mathit{x}}_{\mathrm{a}})\right],\end{array}\end{array}$$

where $\mathbf{K}(=\partial \mathit{y}/\partial \mathit{x})$ represents the Jacobian matrix (n×m), which corresponds to a footprint obtained from the WRF–STILT model, l is the iteration number, and γ is the Levenberg–Marquart parameter fixed at 10 (Chen et al., 2022). The posterior error covariance matrix was calculated with the following equation:

$$\begin{array}{}\text{(7)}& \widehat{\mathbf{S}}=({\mathbf{K}}^T{\mathbf{S}}_{\mathit{\epsilon }}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}{)}^{-\mathrm{1}}.\end{array}$$

The averaging kernel matrix represents the sensitivity of the posterior solution $\widehat{\mathit{x}}$ to the true emission state:

$$\begin{array}{}\text{(8)}& \mathbf{A}={\displaystyle \frac{\partial \widehat{\mathit{x}}}{\partial \mathit{x}}}=\mathbf{I}-\widehat{\mathbf{S}}{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}},\end{array}$$

where I is the identity matrix. The degree of freedom for signal (DOFS), which is the trace of the averaging kernel, indicates the number of independent pieces of information retrieved from the observing system. As a measure of the uncertainty reduction after the inverse analysis, the difference between the prior flux uncertainty and the posterior flux uncertainty relative to the prior flux uncertainty, r, can be used:

$$\begin{array}{}\text{(9)}& r=\left(\mathrm{1}-{\displaystyle \frac{\sqrt{\mathrm{diag}\left[\widehat{\mathbf{S}}\right]}}{\sqrt{\mathrm{diag}\left[{\mathbf{S}}_{\mathrm{a}}\right]}}}\right)\times \mathrm{100}.\end{array}$$

Because we applied weekly and diurnal correction factors from the TIMES model to the hourly footprints in summing them over the STILT run time, we optimized one static emission distribution during the campaign period, assuming that the temporal variation of the emissions followed the TIMES model. Similarly, a single average scaling factor for the large point sources was optimized from the data over the entire campaign period. Considering the numbers of the observation sites and the ΔXCO₂ data available for the inversion (n=654), we aggregated the area source emissions and footprints in the original 1 km × 1 km (0.0083^∘ × 0.0083^∘) grids to $\mathrm{0.025}{}^{\circ }\times \mathrm{0.025}{}^{\circ }$ grids (m=1921). In the sensitivity test, the resolutions of area source emissions and footprints were further lowered to 0.05^∘ × 0.05^∘ grids (m=481) and 0.1^∘ × 0.1^∘ grids (m=121) (Sect. 4.3).

Figure 5 Average CO₂ emission fluxes from (a) ODIAC2020b data in February and March 2016 and (b) MOSAIC data in February and March 2015. (c) The difference between the two datasets aggregated to $\mathrm{0.025}{}^{\circ }\times \mathrm{0.025}{}^{\circ }$ spatial resolution, calculated as (ODIAC − MOSAIC)/(0.5 × (ODIAC + MOSAIC)) × 100. Note that large point sources have been excluded.

To construct the prior error covariance matrix S_a, we compared the ODIAC emission data used as the prior estimate (Fig. 5a) with the MOSAIC emission data (Saito et al., 2022). Although the two databases have similar spatial resolution, it is not the same, and therefore we re-gridded the MOSAIC emission data into the ODIAC grid (Fig. 5b). Then, we aggregated both datasets into the 0.025^∘ × 0.025^∘ grid used for the inverse analysis. Figure 5c shows the difference between the aggregated datasets, calculated as (ODIAC − MOSAIC)/(0.5 × (ODIAC + MOSAIC)) ×100. The difference between the spatial distributions of the ODIAC and MOSAIC data is based solely on area sources because the large point sources are common to the two datasets. Similar large spatial differences also exist among other emission inventories (Gately and Hutyra, 2017) because this kind of emissions database uses geospatial information or physical proxies to allocate the spatial distributions of emissions. Considering the standard deviation of the difference between the two datasets, we set the diagonal elements of S_a to 85 % of the prior emission values. For the scaling factor of the large point source emissions, we set the uncertainty to 15 % based on the temporal variability of monthly liquid natural gas consumed by natural gas-fired power plants of the Tokyo Electric Power Company Holdings, which were available up to March 2016 (https://www.enecho.meti.go.jp/statistics/electric_power/ep002/results_archive.html, last access: 25 July 2023).

The off-diagonal elements of the prior error covariance matrix, which represent the spatial coherence between the prior flux uncertainties in different grid cells, were calculated according to a model of exponential decay with distance between grid cells (e.g., Lauvaux et al., 2016; Lopez-Coto et al., 2020). Thus, the element [i, j] of the prior error covariance matrix was given as

$$\begin{array}{}\text{(10)}& {\mathbf{S}}_{\mathrm{a}}[i,j]={\mathit{\sigma }}_i{\mathit{\sigma }}_j\exp \left(-d_{i,j}/l_{\text{S}}\right),\end{array}$$

where σ_i (σ_j) represents the uncertainty of the emissions in grid cell i (j), d_i,j is the distance between grid cells i and j, and l_S is the spatial correlation length of the prior flux uncertainties. To determine the spatial correlation length of the prior flux uncertainties, we computed semi-variograms of the differences in area source elements between the two emissions datasets in the inversion domain, and then we fitted an exponential model to the semi-variograms with the distance between grid cell pairs limited to 30 km (Mallia et al., 2020). This analysis yielded a correlation length of approximately 10 km (Fig. S4), which is equivalent to that in New York (Pitt et al., 2022) and Salt Lake City (Mallia et al., 2020).

To estimate the measurement error (or model–observation mismatch) covariance matrix S_ε, we used the residual error method of Heald et al. (2004). In this method, the residual errors between the XCO₂ values measured by the EM27/SUN spectrometer and those simulated by WRF–STILT using the prior data were computed, and the variance of the residual over the campaign period was used to represent the diagonal elements of S_ε.

4 Results and discussion

4.1 XCO₂ measurements

During the 2016 Tokyo campaign, the daily minimum XCO₂ values observed by the four spectrometers increased gradually from 403 to 405 ppm as a result of seasonal variation, whereas the daily maximum values showed large day-to-day variation with peaks of up to 415 ppm (Fig. 2a). The XCO₂ values observed at Saitama and Sodegaura from 16 February to 6 April were generally higher than those observed at Tsukuba. To characterize the diurnal variation in XCO₂ at each observation site, we examined the diurnal variation in XCO₂ enhancements (XCO${}_{\mathrm{2}}^{\text{Enh}})$ above the daily XCO₂ baseline. The daily XCO₂ baselines were assumed to be the 5th percentile values of the Tsukuba TCCON measurements throughout each day and to be common to the three sites. We note that the XCO${}_{\mathrm{2}}^{\text{Enh}}$ values were calculated using only the observed XCO₂ values, whereas the ΔXCO₂ values represent the simulations of local XCO₂ enhancement. For days when measurements at Tsukuba were not available (16, 17, 27, and 28 February and 23 March), we used CarbonTracker CT2019B XCO₂ data (see Sect. 3.4). The maximum XCO${}_{\mathrm{2}}^{\text{Enh}}$ value was 9.5 ppm at Saitama and 9.3 ppm at Sodegaura. The average diurnal XCO${}_{\mathrm{2}}^{\text{Enh}}$ value per 15 min bin was calculated for each site using the entire field campaign dataset (Fig. 6). The XCO${}_{\mathrm{2}}^{\text{Enh}}$ values at Tsukuba gradually decreased over time with small standard deviations. This finding reflects the absence of other large emission sources around Tsukuba and the moderate effect of photosynthesis. The measurements at Saitama and Sodegaura showed larger XCO${}_{\mathrm{2}}^{\text{Enh}}$ values than those at Tsukuba, and XCO${}_{\mathrm{2}}^{\text{Enh}}$ values increased over time from approximately 08:00 JST. We note that the high early morning values at Saitama may reflect an airmass-dependent bias. The airmass-dependent variation in XCO₂ is caused by the effects of inaccurate spectroscopic parameters on the retrievals, which vary with the depth of the absorption lines (i.e., airmass) (Wunch et al., 2015). Although this effect is corrected in the GGG2014 software, the error may remain for a large airmass. Although the XCO${}_{\mathrm{2}}^{\text{Enh}}$ variability at Sodegaura was smaller than that at Saitama, some data bins with large standard deviations most likely represent occasional influences of emissions from nearby large point sources.

Figure 6 Average diurnal variations in XCO₂ differences (XCO${}_{\mathrm{2}}^{\text{Enh}})$ from daily background values. These background values were assumed to be common to the three sites and be the 5th percentile value of the Tsukuba TCCON measurements throughout each day. The average XCO${}_{\mathrm{2}}^{\text{Enh}}$ values (open diamonds) and their standard deviations (shading) were calculated for 15 min bins using all data acquired during the campaign period.

When the 2nd (10th) percentile values of the Tsukuba TCCON measurements were used as the daily XCO₂ baseline, the maximum XCO${}_{\mathrm{2}}^{\text{Enh}}$ values were 9.6 (9.4) ppm at Saitama and 9.5 (8.9) ppm at Sodegaura. These changes had little effect on the standard deviations of the mean XCO${}_{\mathrm{2}}^{\text{Enh}}$ values and the pattern of the diurnal variation.

The daily minimum XCH₄ values (Fig. 2b) showed relatively larger temporal fluctuation than those of XCO₂ because of synoptic-scale events (e.g., 22–28 February). Although detailed analysis of XCH₄ data is beyond the scope of this study, in the case study evaluating the WRF–STILT simulation presented in the next section, the XCH₄ values are used.

4.2 XCO₂ simulations

As described in Sect. 3.4, the XCO₂ enhancement (ΔXCO₂) was calculated from the column-averaged footprint and the surface fluxes from area sources, large point sources, and biological activity. Figure S5 in the Supplement shows the ΔXCO₂ values at the three sites separately simulated using area sources, large point sources, and biogenic fluxes. We calculated the contributions of the respective fluxes to the simulated ΔXCO₂ at each site (Table 4). The contribution from area source emissions dominated the simulated ΔXCO₂ for the Saitama and Tsukuba sites. Because the Sodegaura site is located near large point sources, the closest one being approximately 4 km away, the contribution from large point sources was greater at Sodegaura than at Saitama. Because our observations were made from late winter to early spring, the biogenic flux contribution was relatively small, but not negligible, especially at the Saitama site.

Figure 7 Comparison of the XCO₂ observations with the WRF–STILT simulation results (open blue diamonds) at (a) Saitama and (b) Sodegaura. The observations are presented as individual values (open gray circles) and as the 15 min averaged values used for the inversion (open black diamonds).

We compared the XCO₂ data for the forward simulations, which correspond to the XCO₂ simulations from the footprints and the surface CO₂ fluxes based on Eqs. (2)–(4), with the EM27/SUN observations at Saitama and Sodegaura (Figs. 7 and 8). In most cases, the forward simulations captured well the observed temporal variation of XCO₂. However, in some cases they failed to reproduce the diurnal variations. As shown in Fig. 8d, we found it difficult to correctly capture the timing of short-term (<1 h) XCO₂ enhancements, which were probably caused by the plume from large point sources such as the power plants and steel plants located near the Sodegaura site. Similarly, the WRF–STILT simulation for Saitama on 3 March 2016 was not able to capture the XCO₂ enhancement in the late afternoon (Fig. S6a). STILT simulations conducted using ERA5 data and WRF data with different PBL schemes (not shown) showed similar tendencies. Furthermore, even when we changed the emission data from ODIAC to MOSAIC, the discrepancy was not reduced. In addition, we investigated the XCH₄ data, for which a diurnal variation similar to XCO₂ was observed. A WRF–STILT simulation using the Emissions Database for Global Atmospheric Research (EDGAR) version 6 (Crippa et al., 2019) as the CH₄ emission inventory also could not capture the XCH₄ enhancements in the late afternoon (Fig. S6b). However, we cannot rule out the possibility that short-term local sources not included in the prior fluxes may cause the discrepancy between the prior simulations and the observations. Therefore, we attribute this large model–observation discrepancy to errors in the WRF-STILT model, or to the short-term local sources not included in the prior fluxes, or both. These additional simulations indicate that further improvement of the WRF simulation (i.e., assimilation of measurement data) is necessary for more accurate generation of meteorological fields. In our inverse analysis, this large modeling error was considered when setting the measurement error covariance matrix, as described below.

Table 4 Mean fractions of ΔXCO₂ simulated using the three CO₂ fluxes (ΔXCO${}_{\mathrm{2}}^{\text{Area}}$ for area source emission, ΔXCO${}_{\mathrm{2}}^{\text{LPS}}$ for large point source emission, and ΔXCO${}_{\mathrm{2}}^{\text{Bio}}$ for the biogenic flux) to the sum of ΔXCO${}_{\mathrm{2}}^{\text{NPS}}$, ΔXCO${}_{\mathrm{2}}^{\text{LPS}}$, and the absolute value of ΔXCO${}_{\mathrm{2}}^{\text{Bio}}$ for each site.

Figure 8 Comparison of the XCO₂ observations with the WRF–STILT prior (blue) and posterior (red) simulation results for three representative days at (a–c) Saitama and (d–f) Sodegaura. The observations are presented as individual values (open circles) and as the 15 min averaged values used for the inversion (open diamonds).

Figure 9 is a scatter plot between the measured and simulated XCO₂ values over the campaign period; here, the standard deviation of the residual, σ_ε, is 1.31 ppm. As described in Sect. 3.5, the variance of the residual was used as the diagonal elements of S_ε. When the residual between the simulation and observation was more than 3 times σ_ε, the measurements were screened out by greatly increasing the uncertainty. An exponential covariance model in time was selected with a temporal correlation length of 1 h based on the value reported for continuous CO₂ observations in urban areas (Turner et al., 2020).

Table 5 Total CO₂ emissions from the TMA, scaling factors of large point source (LPS) emissions, degree of freedom for signal (DOFS), and prior and posterior XCO₂ differences between the simulations and observations for the different meteorological data, prior emission data, prior uncertainty (σ_a), spatial correlation length of S_a (l_S), temporal correlation length of S_ε (l_t), and spatial resolution of the inversion domain (r_S).

^∗ Data from Sodegaura on 23 March 2016 were excluded.

Figure 9 Scatter plot of observed XCO₂ values and values simulated from the prior (black) and posterior (red) emission fluxes. The mean difference between the simulations and observations (simulation minus observation) with the standard deviation (±1σ) is denoted as δ, and r is the correlation coefficient.

Although the residual error method provides a realistic model–observation mismatch, we also estimated individual uncertainties in our model–observation system, consisting of uncertainties in the measurement data, transport modeling, biogenic flux, and background value. We assumed the uncertainty in measurement data to be the standard deviation of the differences between the EM27/SUN XCO₂ data acquired by side-by-side instruments (Sect. 2). The standard deviation of the bias-corrected XCO₂ differences between the SN38 and SN44 EM27/SUN spectrometers was 0.16 ppm. To estimate the uncertainty in XCO₂ due to the transport modeling error, we ran XCO₂ simulations using the WRF data with different PBL schemes (see Sect. 3.2) and the ERA5 data. The mean biases and the standard deviations of the difference between the EM27/SUN measurements and the STILT simulations from the prior fluxes are listed in Table 5 (Prior XCO₂ difference). Whereas there was no large difference in the standard deviation among the three simulations using the WRF data (1.31–1.40 ppm), the standard deviation for the simulation using ERA5 was 2.74 ppm, more than 1 ppm larger than that of the simulations using WRF. This large value was because the ΔXCO₂ values at Sodegaura on 23 March 2016 simulated from the ERA5 data showed a rather large peak (∼20 ppm) caused by incidental contamination from the nearby large point sources that was not present in the actual measurement data. When the data for that site and day were excluded, the standard deviation decreased to 1.85 ppm. The XCO₂ uncertainty resulting from transport modeling, estimated as the standard deviation of the differences between the XCO₂ values simulated using the WRF and ERA5 meteorological fields, was 1.65 ppm. To estimate the uncertainty in XCO₂ resulting from the biogenic flux error, we calculated XCO₂ values over the campaign period for four types of biogenic fluxes (VISITc and three others) with differing spatial and temporal resolutions (Table S1 in the Supplement) but with other input parameters unchanged. The Simple Biosphere Model version 4.2 (SiB4, Haynes et al., 2021) and the Biosphere model integrating Eco-physiological And Mechanistic approaches using Satellite data (BEAMS, Sasai et al., 2005) are both terrestrial biosphere models, whereas CarbonTracker version CT2019B (Jacobson et al., 2020) is from a data assimilation system in which the biogenic and oceanic fluxes are optimized. The average standard deviation across the XCO₂ values calculated using the four biogenic fluxes, 0.09 ppm, was regarded as the XCO₂ uncertainty resulting from the biogenic flux. When Tsukuba TCCON data were not available and CarbonTracker data were used instead, the uncertainty in the background value rose. We assumed that the XCO₂ uncertainty resulting from the background value was represented by the standard deviation of the XCO₂ difference between the Tsukuba TCCON data and the CarbonTracker data and estimated the uncertainty as 0.72 ppm (see Sect. 4.1). These evaluations revealed that the uncertainty in transport modeling was dominant, followed by the uncertainty in background value.

Figure 10 shows the mean ΔXCO₂ contribution from area source emissions during the field campaign, which was calculated from the absolute value of the ΔXCO₂ difference between the urban and Tsukuba sites. We limit the domain for estimating area source emissions to the urban domain indicated by the magenta rectangle in Fig. 10 (hereinafter referred to as the “inversion domain”), where the contributions of each grid cell to the modeled ΔXCO₂ are relatively large.

4.3 Posterior fluxes

Figure 11a shows the posterior area source CO₂ emissions, estimated using the settings for the reference inversion (i.e., case #0 in Table 5 using the ODIAC as the prior data and meteorological fields from the WRF model with the MYJ PBL scheme for footprint calculations). The total DOFS from the reference inversion was 6.49, of which 5.73 is for spatially resolved emissions and 0.76 is for the large point source emissions. The spatial pattern of the optimized emissions still largely resembled the prior estimate pattern (Fig. 11a). In large parts of Tokyo and Kanagawa, the emissions were revised downward, whereas in Saitama, Ibaraki, and northern Chiba, the emissions became larger. Because the mean bias in XCO₂ values simulated from the prior emission flux was originally small, the emissions from the central TMA region became smaller than the prior values, and the emissions from the other regions became larger than the prior values. The spatial distribution of the changes from the prior flux (Fig. 11b) was partly in agreement with the spatial differences between the MOSAIC and ODIAC emission data (Fig. 5c). Because the locations of large point sources were corrected in the prior emissions (Sect. 3.3), the difference in the spatial distribution between the prior and posterior emissions may be due to the non-representativeness of the spatial proxy (i.e., night lights) used in the ODIAC data. As an indication of the efficiency of the inversion, we evaluated to what extent the differences between the XCO₂ simulations and observations were improved by using the posterior fluxes. The XCO₂ values simulated from the posterior fluxes were in better agreement with the observations than those simulated from the prior fluxes (Fig. 8). The mean bias in XCO₂ simulations against observations decreased from 0.30 to −0.03 ppm, but the RMSE decreased only slightly, from 1.31 to 1.21 ppm (Fig. 9). This slight RMSE reduction is because the emission distribution was estimated on a monthly basis, whereas the individual model–observation discrepancies were governed by the transport modeling error. Next, we compared the estimated total emissions in the TMA with the emission inventories. The total emissions correspond to the domain-aggregated emission flux during the campaign period (i.e., from February to March 2016). Figure 12 shows the total emissions calculated from the prior flux and the posterior flux in the reference inversion. The error bars (uncertainties at the 95 % confidence level) of prior and posterior total emissions are based on the respective error covariance matrices and were obtained by summing the emission uncertainties in each grid cell and the uncertainty of the large point source emission in quadrature. The posterior large point source emissions were adjusted downward by 14.4 % compared with the prior emissions (i.e., scaling factor of 0.856), and the posterior area source emissions were adjusted upward by 10.4 %. Consequently, the difference between the prior and posterior total emissions was approximately 1 %. Although the change in the total emissions was relatively small, the inversion led to a reduction of the uncertainty in the total emissions by a factor of ∼2 (i.e., the uncertainty at the 95 % confidence level decreased from 11.3 % to 5.2 %).

Figure 10 Mean contribution of each grid cell to the ΔXCO₂ values simulated from prior area source emissions over the campaign period. The CO₂ emissions were optimized for the domain within the magenta rectangle.

Figure 11 (a) Area source CO₂ emission fluxes in the TMA in the reference inversion (case #0 in Table 5) combining posterior (within the magenta rectangle) with prior (outside the rectangle) CO₂ emissions and (b) the difference between the posterior and prior emissions. The dotted lines show the administrative boundaries.

Here, we present the results of emission estimates obtained for cases with different inversion settings (Table 5): case #1, large point source emissions fixed; cases #2a–c, footprints calculated from different meteorological fields used; cases #3a and #3b, prior uncertainty halved or doubled; cases #4a and #4b, spatial correlation length of S_a changed; cases #5a and #5b, temporal correlation length of S_ε changed; case #6, EDGAR version 6 (0.1^∘ × 0.1^∘ spatial resolution) without large point source correction used as the prior estimate (Fig. S7); and cases #7a and #7b, spatial resolution of the inversion domain coarsened to 0.05^∘ or 0.1^∘ (i.e., 2 or 4 times the reference case). For the case #6 and #7 inversions, the prior uncertainty and the spatial correlation length were re-determined as described earlier. The total emissions, scaling factor of large point source emissions, and ΔXCO₂ bias between the simulation and observations and its standard deviation for each case are summarized in Table 5. For case #1, the inversion in which the large point source emissions were fixed, both the mean bias and the standard deviation of the posterior XCO₂ simulations against observations were equivalent to those of the reference inversion (case #0). Although the scaling factor of large point source emissions for the reference inversion was 0.856, total emissions in case #1 were 5.3 % larger than those in case #0. The posterior XCO₂ simulation results obtained with different meteorological fields (cases #2a–c) indicated that the biases and standard deviations were improved compared to the prior XCO₂ simulations, irrespective of the meteorological field. Among them, use of the WRF model with the MYJ scheme resulted in the smallest standard deviations for not only the prior but also the posterior XCO₂ simulations. When the prior uncertainty (cases #3a and #3b), its correlation length (cases #4a and #4b), and the temporal correlation length of S_ε (cases #5a and #5b) were changed, the mean biases and the standard deviations of the posterior XCO₂ simulations were comparable to the reference inversion. However, the prior uncertainty had a larger impact on the total emission estimates than the spatial correlation length. The inversion using EDGAR as the prior emission inventory (case #6) resulted in a posterior XCO₂ simulation with a low bias of 0.16 ppm and a standard deviation of 1.27 ppm; this simulation underestimated the total emissions by 15.8 % compared with the reference inversion. This result implies that the use of emission data with a low spatial resolution introduces additional uncertainty into XCO₂ modeling. Similarly, reducing slightly the spatial resolution of the ODIAC data (cases #7a and #7b) degraded both the mean bias and the standard deviation of the posterior XCO₂ simulations. In the case of the reference inversion, the number of measurement data points was considerably smaller than the number of grid cells whose emissions were optimized. Although the number of grid cells with a spatial resolution of 0.05 and 0.1^∘ was equivalent to or lower than the number of measurement data points, respectively, the total DOFS slightly decreased (to 5.84 for 0.05^∘ and 5.05 for 0.1^∘). This was due to the changes in the prior uncertainty and the spatial correlation length. The posterior total emissions did not differ greatly from those of the reference inversion. Figure 12 shows the ensemble mean of the total emissions and its uncertainty at the 95 % confidence level, estimated from the scatter of these inversion results. For comparison, the total emissions from the original ODIAC data and the original and LPS-corrected MOSAIC data are also displayed. The ensemble mean total emissions are in agreement with the original and LPS-corrected MOSAIC emissions within the uncertainty of the ensemble inversions.

Figure 12 Total CO₂ emissions in the TMA calculated from the posterior emission fluxes (red), ODIAC data (blue), and MOSAIC data (green). The posterior emission fluxes are shown for the reference inversion (case #0 in Table 5) and the ensemble mean of all cases listed in Table 5. The error bars for the prior and posterior emission fluxes are the respective estimated uncertainties, whereas the error bar for the ensemble mean is the standard deviation. For the ODIAC and MOSAIC data, both original and LPS-corrected total emissions are shown.

We compared our results with those of a previous CO₂ inversion study for the TMA (Pisso et al., 2019; Babenhauserheide et al., 2020) and with annual emissions in fiscal year (FY) 2015 (April 2015 to March 2016) reported by each administrative division in the Tokyo Metropolis. Here, we calculated the total CO₂ emissions in the Tokyo Metropolis by integrating emissions in the grid cells within its administrative boundaries. Additionally, because our inversion domain included almost the whole area of the Tokyo Metropolis and of Kanagawa, Chiba, and Saitama prefectures, the total emissions from these four administrative divisions (referred to as “southern Kanto”) were also calculated. The total emissions estimated by our reference inversion were 56.6 Mt-CO₂ yr⁻¹ for the Tokyo Metropolis and 277.8 Mt-CO₂ yr⁻¹ for southern Kanto (Fig. S8), and these emissions are smaller by 29 % and 50 %, respectively, than those estimated by Pisso et al. (2019). We note that although Pisso et al. (2019) estimated mean emissions for 2005–2009, the difference between the FY2015 emissions and the FY2005–2009 mean emissions reported by the Tokyo Metropolis is less than 1 % (https://www.kankyo.metro.tokyo.lg.jp/en/climate/index.files/Tokyo_GHG_2019.pdf, Pisso et al., 2019) and this study uses comparable Lagrangian transport models to calculate atmospheric transport; however, there are several differences, including the type of observational data (in situ vs. column), the prior emission fluxes (EDGAR vs. ODIAC), the meteorological fields for driving the transport model (ERA-Interim vs. WRF based on GPV-MSM), and the spatial resolution of emission estimates (20 km × 20 km vs. 3 km × 3 km). Our sensitivity analysis shows that changing the prior fluxes, meteorological field, and emission estimation resolution to roughly match Pisso et al. (2019) did not produce a result substantially different from the emission estimation result of the reference inversion. We thus concluded that the improved accuracy of emission estimates in our study may be due to the use of columns as observational data. Column data are less susceptible to the effect of PBL height changes that are difficult to simulate in transport models and have information on a larger area of emissions due to the difference in wind direction at each altitude. Babenhauserheide et al. (2020) estimated CO₂ emissions of 256 ± 77 Mt-CO₂ yr⁻¹ for the urban area around Tokyo. Our emission estimate for southern Kanto was in reasonable agreement with the result of Babenhauserheide et al. (2020), although the comparison is not exact because of the discrepancy in the areas where the CO₂ emissions were calculated. The total emissions in FY2015 reported from each administrative division were 60.3 Mt-CO₂ yr⁻¹ for the Tokyo Metropolis and 250.6 Mt-CO₂ yr⁻¹ for southern Kanto (Table S2); these values show remarkable agreement with our posterior estimates from the reference inversion. Furthermore, our posterior estimate for the Tokyo Metropolis lies between the ODIAC and MOSAIC inventory data. The relationship between our posterior estimate and the inventory data for southern Kanto is similar to that for the TMA shown in Fig. 12, because southern Kanto includes most of the TMA as defined in this study. Thus, these comparisons demonstrate that our top–down approach was able to properly constrain CO₂ emissions in this urban area.

5 Conclusions

We conducted a field campaign to estimate CO₂ emissions in the TMA from February to April 2016 with two EM27/SUN spectrometers deployed at sites in Saitama and Sodegaura and the Tsukuba TCCON spectrometer. The XCO₂ values at Saitama and Sodegaura exhibited large enhancements compared with those at Tsukuba, and the mean diurnal variation of the enhancements showed a tendency to increase toward evening. The Lagrangian transport model STILT, which was driven by WRF meteorological fields generated at a horizontal resolution of ∼1 km, was used for simulating the XCO₂ enhancements resulting from anthropogenic (area source and large point source) emissions and biogenic fluxes. As with the prior fluxes, the anthropogenic emissions from the ODIAC dataset were corrected by replacing the locations and emission magnitudes of large point sources with inventory data, whereas the biogenic flux from VISITc was downscaled using GVF data. We found that, for the TMA, the WRF model with the MYJ PBL scheme and the RUC land surface model yielded optimal results with regard to both wind fields and the XCO₂ simulations. The XCO₂ forward simulation results using the prior fluxes highlight several factors that should be considered when designing an observation campaign or an operational network for ground-based column measurements for estimating urban emissions. Although the XCO₂ forward simulations generally showed good agreement with the observations, the comparison between the simulations and observations demonstrated some limitations in the modeling capability. As described in Sect. 4.2, in some cases the simulations failed to reproduce the diurnal variation and to capture the plume from nearby large point sources, possibly because of the transport modeling error or the short-term local sources not included in the prior fluxes (Figs. 8d and S6). Assimilating meteorological measurement data such as AMeDAS into the WRF calculation would be one way to reduce the modeling error. Additional wind lidar observations would be useful to better constrain wind fields and PBL as a whole (Deng et al., 2017). However, it is a great challenge to simulate local plumes from large point sources. In a previous study, we conducted simultaneous measurements of XCO₂ and wind data with the EM27/SUN instruments and a Doppler lidar, respectively, co-located close to a thermal power plant in Japan (Ohyama et al., 2021). Because not even the simulation using the measured wind data and a simple dispersion model could reproduce the timing of the observed XCO₂ enhancement, we decided to adjust the wind directions as part of the optimization of emission fluxes. At the Sodegaura site, where there are two large point sources within 10 km and two more within 15 km, the contribution from the large point sources in the TMA to the simulated ΔXCO₂ is equivalent to the contribution from area sources (Table 4 and Fig. S4). These findings suggest that, for the purpose of estimating emissions from the entire city, the locations of the EM27/SUN instruments should be selected to avoid proximity to large point sources or, through consideration of the dominant wind direction, to minimize the influences from large point sources.

Using these observational and modeling approaches along with their uncertainties, we developed an urban area-scale inversion system to estimate spatially resolved CO₂ emission at >3 km resolution and a suitable scaling factor for large point source emissions. The posterior CO₂ flux reduced both the mean bias and the standard deviation of the differences between the XCO₂ simulations and observations. Whereas the posterior total CO₂ emissions in the TMA from the reference inversion were consistent with those from the prior estimate with ∼1 %, the posterior uncertainty was halved compared with the prior uncertainty. The ensemble mean of the posterior total CO₂ emissions agreed with the LPS-corrected ODIAC (prior) and MOSAIC data within the posterior uncertainty at the 95 % confidence level estimated from the ensemble scatter. We conclude that the EM27/SUN data could constrain urban-level CO₂ emissions and partially resolve the spatial distribution at monthly scale. Because few EM27/SUN instruments were available for the 2016 Tokyo campaign, we deployed only two EM27/SUN instruments with consideration of the prevailing wind direction. The actual wind direction varied more than expected, with the result that about 1 month of data showed a wide range of sensitivity, as shown in Fig. 10. The deployment of additional instruments would increase our sensitivity to emissions and thus the DOFS. This would enable more frequent (i.e., bimonthly or weekly) emission estimates. In addition, more instrument locations would also help to constrain the background. We plan to construct operational observation sites with EM27/SUN spectrometers in central Tokyo and the TMA suburbs. These data not only will facilitate the operational estimation of CO₂ emissions in the TMA, thereby helping to verify emission reduction efforts, but also will validate GHG data from future satellite missions with small footprints and a wide swath width, such as Japan's GOSAT-GW (Global Observing SATellite for Greenhouse gases and Water cycle; https://gosat-gw.nies.go.jp/en/, last access: 25 July 2023) and ESA's CO2M (Copernicus Anthropogenic Carbon Dioxide Monitoring mission; https://www.esa.int/ESA_Multimedia/Images/2022/03/CO2M, last access: 25 July 2023).