Estimation of CO₂ fluxes in the cities of Zurich and Paris using the ICON-ART CTDAS inverse modelling framework

Observation-based estimation of urban CO₂ emissions can help cities track their pathway to net zero emissions, a goal many cities worldwide have adopted. While mesoscale atmospheric transport models are an effective component in inversion systems estimating country-level emissions, their use in urban-scale inversions presents a significant challenge. Here, we present one-year flux inversion results with the mesoscale ICON-ART atmospheric transport model for two cities with contrasting size and topographic complexity: Zurich and Paris. Inversions were performed with an ensemble square root filter, assimilating observations from a dense rooftop CO₂ sensor network in Zurich and from a tall tower network in Paris. The inversion framework optimized gridded anthropogenic and biospheric fluxes, along with background mole fractions from eight inflow regions. Prior anthropogenic emissions were based on detailed inventories provided by local authorities. In Zurich, the inversion resulted in a posterior annual anthropogenic emission of 1012.3±38.8 kt yr⁻¹, representing approximately a 30 % reduction compared to the prior, with the most significant decreases during winter periods of elevated ambient temperatures. In contrast, the posterior fluxes in Paris remained close to the prior, with an annual emission of 3580.0±101.9 kt yr⁻¹, which is 7 % higher than the prior. This comparison highlights the influence of city-specific factors – such as topography, city size, and observational network – on the inversion system performance. Furthermore, our findings demonstrate the potential of mesoscale models to refine urban emission estimates, offering valuable insights for policymakers and researchers working to improve emission inventories and advance urban climate strategies.

1 Introduction

Inversion of CO₂ fluxes using atmospheric transport models is a well-established approach that was originally applied at global scale to constrain the global carbon budget (Gurney et al., 2002; Davis et al., 2021; Broquet et al., 2013; Monteil et al., 2020; van der Woude et al., 2023). With the development of regional measurement networks and advances in high-resolution modeling, this approach has since been extended to continental and national scales. However, despite growing interest, only a limited number of inverse modeling studies have been performed at the urban scale, for example (Lauvaux et al., 2016; Kunik et al., 2019; Sargent et al., 2018; Nalini et al., 2022). Urban inversions are gaining increasing attention, as cities are major contributors to anthropogenic emissions (44 %, Seto et al., 2014), making accurate emission estimates at this scale essential for supporting climate action plans and verifying reported emission reductions.

In this study, we present results from year-long CO₂ flux inversions conducted for two European cities, Zurich and Paris. This work is part of the European ICOS Cities project, which aimed to develop and evaluate different CO₂ emission monitoring systems in three pilot cities of contrasting size and topographic complexity – Paris, Munich and Zurich. The primary objective of our study is to generate robust estimates of anthropogenic CO₂ emissions in Zurich and Paris, using high-resolution atmospheric transport simulations in combination with dense urban CO₂ observation networks.

Previous studies estimating urban CO₂ emissions have highlighted several key challenges. One major difficulty is to differentiate between CO₂ enhancements caused by local sources and those resulting from inflow from surrounding regions (Lauvaux et al., 2016; Sargent et al., 2018; Broquet et al., 2013). This is because the enhancements in CO₂ mole fractions due to biospheric and anthropogenic sources within the domain are often of similar amplitude as the variations in the background levels. As a result, it is crucial to include observation sites located outside – and ideally upwind – of the urban area to better constrain background conditions. In this study, we address this issue by jointly optimizing background mole fractions along with anthropogenic and biospheric fluxes. Without this simultaneous optimization, any bias in the background would propagate directly into the estimated fluxes, compromising the accuracy of the inversion.

Another challenge is the separation of anthropogenic and biospheric contributions (Lauvaux et al., 2016; Sargent et al., 2018; Kunik et al., 2019; Wu et al., 2018), as their atmospheric signals often overlap in both space and time. Different attempts have been made in previous studies to separate biospheric fluxes from anthropogenic emissions. Additional challenges arise at the urban scale. For example, Lian et al. (2023) showed for Paris that coarse land-use vegetation data can underestimate urban biosphere activity considerably. They also mention the lack of eddy covariance measurements over urban vegetation to validate prior biospheric flux estimates. Sargent et al. (2018) demonstrated for Boston that the urban biosphere can take up more than half of the anthropogenic CO₂ signal during summer afternoons, demonstrating the importance of constraining biospheric fluxes alongside anthropogenic emissions. Similarly, synthetic inversions over Salt Lake City and Indianapolis emphasized that prior assumptions on biospheric fluxes strongly affected posterior emissions and that dense, strategically placed measurement networks are crucial to disentangle overlapping signals (Kunik et al., 2019; Lauvaux et al., 2016).

A final challenge relates to the prior information on fluxes and their uncertainties (Broquet et al., 2013; Lauvaux et al., 2016; Sargent et al., 2018; Kunik et al., 2019). In our case, we could benefit from very detailed city inventories provided by the authorities of Zurich and Paris.

The measurement networks used in the inversions differed significantly between Zurich and Paris, reflecting their contrasting geographic and urban characteristics. In Paris, located in the flat Île-de-France region, a tower-based network of 9 sites with high-precision instruments was operated. The network was designed to measure both upwind background mole fractions and downwind increments due to urban emissions. Depending on prevailing wind direction, several sites alternated between upwind and downwind roles.

In contrast, Zurich's complex topography – characterized by intersecting valleys and surrounding ridges – makes it difficult to consistently define upwind and downwind locations. To address this, a denser network of 13 mid-cost rooftop sensors was installed within the city to monitor urban CO₂ enhancements, complemented by three high-precision instruments mounted on towers outside the urban area to provide background measurements.

For the inversion, we used the mesoscale atmospheric transport model ICON-ART (Zängl et al., 2015; Rieger et al., 2015; Schröter et al., 2018) coupled with the CarbonTracker Data Assimilation Shell (CTDAS) (Peters et al., 2005; van der Laan-Luijkx et al., 2017). This coupling of ICON-ART and CTDAS, originally developed for methane inversions at the European scale (Steiner et al., 2024b), was extended in this study to run CO₂ inversions over urban domains and to modify biospheric fluxes “online” (i.e., during the simulation, Jähn et al., 2020). The ICON-ART simulations were performed at a spatial resolution fine enough to resolve the main topographic features of the Zurich area, including the Limmattal and Glattal valleys and the surrounding ridges, which rise 100–400 m above the valley floors.

In this study, we adress the aforementioned challenges by applying a high-resolution inversion framework that jointly optimizes anthropogenic emissions, biospheric CO₂ fluxes, and background concentrations for two cities with markedly different characteristics. The comparison spans a relatively straightforward case – Paris, a large, isolated city in flat terrain – and a more complex scenario – Zurich, a mid-sized city embedded in mountainous terrain and surrounded by other urban agglomerations. By comparing results from Zurich and Paris, we explore how sensor network design, atmospheric transport, and prior flux uncertainties influence inversion performance.

2 Data and Methods

2.1 Observations

Atmospheric CO₂ dry air mole fraction measurements for both cities were obtained from ground-based station networks (see Fig. 1). For brevity, we will refer to dry air mole fractions simply as mole fractions. The measurement data were first aggregated to hourly values and then averaged over afternoon hours (11:00–16:00 UTC, i.e., 12:00–17:00 local time) before being used for model evaluation and flux inversions. It is common practice to consider only daytime or afternoon measurements when comparing simulated and observed mole fractions, as stable nocturnal boundary layers are challenging to simulate and often lead to large model-observation mismatches (Gerbig et al., 2003).

Figure 1 Locations of the measurement stations comprising the Zurich (a) and Paris (b) CO₂ observation networks. Tower sites in both cities are labeled with their acronyms. Photo of the Zurich Schule Milchbuck station by Pascal Rubli; photo of the Paris Saclay station by Michel Ramonet. Basemap tiles for panels  (a, b) © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

In Zurich, the observation network comprised 13 rooftop measurement sites inside the city equipped with mid-cost sensors. Here we use the term mid-cost sensors to distinguish them from high-precision analyzers (e.g., cavity ring-down spectroscopy) and from low-cost sensors such as those deployed in the Carbosense network (Müller et al., 2020). As described in (Grange et al., 2025), three different models of non-dispersive infrared (NDIR) sensors were deployed in Zurich's mid-cost sensor network called ZiCOS-M. Different from the sensors in Zurich's low-cost sensor network ZiCOS-L, which was operated in parallel (Creman et al., 2025), the mid-cost sensors had higher sensitivity, were mostly operated in temperature-controlled rooms, and were calibrated daily by supplying calibration gas from two reference gas cylinders. Their accuracy was about 1 ppm (Grange et al., 2025), which is one order of magnitude better than the accuracy of the low-cost units. Three background sites were located outside the city. Two of the background sites were equipped with high-precision instruments, one with a mid-cost sensor. An overview of the sites is provided in Table 1. Hardau II (hard) is a central site with a 20 m long mast mounted on top of a 95.3 m building. The sampling at this site occurs at a much higher altitude compared to the other rooftops in the network (Grange et al., 2025). This site, which features additional monitoring activities including an eddy covariance system, is highlighted in Fig. 1 as high-rise to indicate its elevated inlet height. The mid-cost sensor network consisted of 19 sites inside the city, but six of these were discarded because they either measured at street-level or were influenced by local sources. Both situations cannot be captured adequately by a mesoscale atmospheric transport model. Data from the stations were subject to basic quality control before being used in the inversion. A full description of the Zurich sensor network, calibration strategy and data processing is provided in Grange et al. (2025). The inversion period spanned a full year from September 2022–August 2023, preceded by a two-week spin-up in August 2022.

Table 1 ICOS Cities sites in Zurich with mid-cost sensors installed on rooftops and 2 background sites (Beromünster and Laegern–Hochwacht) equipped with high-precision instruments. Elevation refers to height above sea level, and inlet height to height above ground level.

In Paris, all measurements used in this work were recorded at nine tower sites equipped with high-precision instruments (Doc et al., 2024). Most of the sites were located outside the city, with three measuring upwind mole fractions and three measuring downwind ones, while the remaining three were located within the city. The sites are listed in Table 2. The inversion for Paris also covered a full year, but for the period between January–December 2023. A detailed description of the measurement network and data handling procedures is available in (Doc et al., 2024). A mid-cost CO₂ sensor network was recently also established in Paris (Lian et al., 2024). These sites were excluded from this study due to upgrades and expansion under the ICOS Cities project, but they will be available for future studies.

Table 2 Sites from the tower network in Paris. All sites are equipped with high-precision instruments. Elevation refers to height above sea level, and inlet height to height above ground level.

The dense rooftop sensor network in Zurich offers detailed spatial coverage and high sensitivity to emission sources within the city. In contrast, the tall tower network in Paris is only sensitive to emissions integrated over larger portions of the city.

2.2 ICON-ART model

ICON-ART (ICOsahedral Nonhydrostatic model with Aerosols and Reactive Trace gases) is a mesoscale meteorology and atmospheric transport model that consists of two components: the climate and weather prediction model ICON (Zängl et al., 2015), developed by the German Weather Service (DWD) and the Max Planck Institute for Meteorology, and ART (Rieger et al., 2015; Schröter et al., 2018; Hoshyaripour et al., 2025), mainly created by the Karlsruhe Institute of Technology, for simulations of passive and chemically reactive tracers. ICON is a versatile model that can be run from global to regional and even to sub-kilometer scale. It operates on a semi-structured grid with triangular grid cells and offers options for online (through regional grid refinement) and offline nesting. Here, we use ICON in limited-area configurations with offline nesting. ICON-ART is a fully coupled model jointly simulating weather and atmospheric tracer transport in a consistent way (Baklanov et al., 2014). ICON and ART have recently been released under a permissive open source license and are now maintained and developed by a broader consortium of German and Swiss research partners and weather services.

Model setup

We ran separate CO₂ simulations for the two cities, each nested offline within a larger central European domain. The central European domain was simulated at 6.5 km resolution using grid R3B8 (see Fig. 2), while the inner domains were run at a resolution of 0.5 km for Zurich (R19B9 grid) and 1 km for Paris (R5B10 grid). Zurich required a higher resolution due to its complex topography and smaller size. The European simulation provided the initial and boundary conditions for the two nested domains and was itself nested into the Copernicus Atmospheric Monitoring Service (CAMS) global inversion-optimized CO₂ simulation v24r3, which is based on assimilation of satellite observations (Chevallier et al., 2010; Copernicus Atmosphere Monitoring Service (CAMS), 2020). Meteorological initial and boundary conditions for the European simulation were taken from the ERA5 global reanalysis of the European Centre for Medium Range Weather Forecasts (ECWMF) (Hersbach et al., 2020). The intermediate European simulation was necessary because of the significant spatial resolution difference between the CAMS and ERA5 global products and the high-resolution simulations for the two cities. To keep the simulated meteorology close to the analyzed meteorology, the ICON-ART meteorological fields were weakly nudged towards the ERA5 data as described in Steiner et al. (2024b).

Figure 2 Model domains used in this study. The outermost domain covers Central Europe at 6.5 km resolution. The insets show the two nested domains centered on the Zurich and Paris metropolitan areas simulated at a resolution of 0.5 and 1 km, respectively. The topography is shown as background shading. The black contour denotes the limits of the city of Zurich. For Paris, the black contour corresponds to the Île-de-France, a region covering the agglomerations of Paris. The city limits are represented by the green contour.

The same parameters were also used in the study of Stagakis et al. (2025), which compared four different biospheric flux models applied over the city of Zurich with each other and against observations of respiration fluxes, leaf area density and sap flow in urban parks. This limited and partially indirect evaluation (e.g. sap flow used as a proxy of GPP) showed VPRM to perform equally well as other, more complex biospheric models. A more extensive evaluation of different variants of VPRM including urban VPRM (Hardiman et al., 2017) against measurements in Zurich and Munich is currently ongoing. It is clear that further developments are needed to improve biospheric flux models for urban areas, since they have mostly been developed for, and tuned to, natural environments.

All simulations were run with 60 vertical layers. The time step was 50 s in the European domain, 10 s for Paris, and 5 s for Zurich. Model output was saved hourly. The Tiedtke–Bechtold convection scheme (inwp_convection = 1) was used throughout, but in the city domains, only shallow convection was enabled (shallowconv_only = .TRUE.), while in the European domain both shallow and deep convection were active. The configuration for the city domains closely followed the setup of the Swiss weather service for its operational forecasts at 1 km resolution.

Assimilation of CO₂ observations was only performed in the two nested domains, not in the European domain. Any potential CO₂ mole fraction biases in the European run, which served as boundary conditions for the regional domains, were corrected for by the inversion as described later. To estimate the magnitude of these biases, the simulated mole fractions were compared against measurements from the European Integrated Carbon Observation System (ICOS) (Yver-Kwok et al., 2021). The results will be presented in Sect. 3.1.1.

Different gridded emission inventories were used depending on the simulation domain. For the European domain, the TNO-GHGco inventory (Super et al., 2020) for the year 2021 at about 5 km resolution was used. This inventory relied on the 2023 official reporting of 2021 emissions from the AVENGERS project, except for shipping emissions, which were sourced from the 2021 TNO-GHGco inventory used in the CoCO₂ project, reflecting an earlier reporting version for that sector. For the Zurich domain, three inventories were combined: the TNO-GHGco inventory was used for regions outside Switzerland; within Switzerland, the 2020 Swiss national inventory at 100 m resolution was applied; and for the city of Zurich specifically, this was replaced by a more detailed 2020 inventory as described in Brunner et al. (2025). The city inventory, originally provided as point, line, and area sources, was rasterized onto the ICON model grid before merging with the other inventories (see Fig. 3a). The three inventories were mapped to the ICON grid and merged into a single dataset using the Python package emiproc (Lionel et al., 2025). For the Paris domain, the TNO-GHGco inventory was merged with a 500 m resolution inventory for the Île-de-France (see Fig. 3b) provided by the regional air quality agency AIRPARIF for 2022.

Figure 3 Annual mean anthropogenic CO₂ emissions in the Zurich (a) and Paris (b) domains in kt yr⁻¹. The Zurich data are based on the Swiss national and Zurich inventories, along with TNO-GHGco inventory covering bordering regions of Germany. The Paris data are from AIRPARIF inventory for Île-de-France region, with the remaining area covered by TNO-GHGco inventory. City borders are shown with black lines, and lakes in the Zurich domain are outlined in white.

The temporal emission profiles were not always available, necessitating the use of generic time profiles in some cases. To account for these limitations and avoid introducing biases into the optimization, relatively large uncertainties were prescribed for all anthropogenic and biospheric sources, allowing the system to adjust the fluxes based on observations. At the same time, spatial correlations were introduced in the prior error covariance matrices to limit the degrees of freedom and prevent overfitting.

Biospheric fluxes were computed online using the Vegetation Photosynthesis and Respiration Model (VPRM, Mahadevan et al., 2008) integrated into ICON-ART. As input, it requires shortwave radiation and two-meter temperature, which were provided by ICON, as well as satellite measurements of two indices, the enhanced vegetation index (EVI) and the land surface water index (LSWI) obtained from the MODIS instrument (Vermote, 2015). These indices describe the influence of phenology and water stress on photosynthetic uptake of CO₂ by vegetation. VPRM independently simulates gross photosynthetic production (GPP) and ecosystem respiration (RE) for seven vegetation classes. The net ecosystem exchange (NEE) is then the difference RE–GPP. The vegetation class-specific parameters of VPRM were taken from Table F1 of Glauch et al. (2025), obtained through an optimization procedure that compared VPRM simulated fluxes with Eddy covariance flux tower observations in Europe.

Another important input for calculating biospheric fluxes is land cover, which provides information about the spatial distribution and type of vegetation in the model domain. We used the CORINE 2018 land cover dataset from the European Environment Agency's Copernicus Land Monitoring Service at 100 m resolution (European Environment Agency (EEA), 2018). However, CORINE fails to resolve small vegetation patches or individual trees in urban areas; thus, a high-resolution vegetation dataset was created for Zurich as described in Brunner et al. (2025) and merged with CORINE data outside the city. Since no such tailored high-resolution product was developed for Paris, vegetation cover in Paris is likely underestimated. A 10 m resolution version of the land cover map for Zurich is presented in the Supplement.

Land cover in ICON-ART is represented using a tile approach, where multiple land cover types are included in a grid cell proportional to their fractional area. We used the default setting of three tiles, i.e., only the three dominating land cover types were represented. In metropolitan areas, these typically consist of one urban land cover (non-vegetated) and two of the seven vegetation classes considered by VPRM, often grassland and deciduous forest.

2.3 Inversion approach

2.3.1 CTDAS inversion framework

To estimate CO₂ fluxes, we used the Carbon Tracker Data Assimilation Shell (CTDAS), which is a flexible framework that employs an Ensemble Square Root Filter (EnSRF) to optimize a state vector of flux scaling factors and other elements (Peters et al., 2005; van der Laan-Luijkx et al., 2017). CTDAS has been coupled with various Eulerian and Lagrangian atmospheric transport models, including ICON-ART (Steiner et al., 2024b). The ensemble approach requires the model to simulate a large ensemble of CO₂ tracers, each representing a different perturbation of the state vector elements being optimized.

Our setup closely followed that in Steiner et al. (2024b). In short, we jointly optimized scaling factors for fluxes and background CO₂ mole fractions over individual 1 week assimilation windows. Each week was optimized twice: first with observations from the current week and a second time with observations from the following week. This assumes that CO₂ emitted in the current week also influenced CO₂ in the next week.

Estimating the fluxes for a full year requires 52 individual assimilation cycles, each cycle consisting of three week-long simulations (except for the first cycle, which simulates only two weeks). A schematic illustrating the workflow and information transfer between weeks and assimilation cycles was presented in Steiner et al. (2024b). Figure 4 presents an alternative view illustrating the ensemble of CO₂ mole fractions simulated by the system at an arbitrary observation location.

Figure 4 Schematic representation of CTDAS assimilation cycles (a) 1,  (b) 2 and  (c) 3. Black dots indicate hourly daytime (11:00–16:00 UTC) observations which were averaged daily and assimilated during each cycle, while gray dots show all other observations. Each panel shows observations only for days assimilated in the corresponding cycle. λ_i denotes the state vector parameters (flux scaling factors) for cycle i, ${\mathit{\lambda }}_i^{a\mathrm{1}}$ and ${\mathit{\lambda }}_i^{a\mathrm{2}}$ indicate once and twice optimized scaling factors, respectively, following the notations in Peters et al. (2005) and Steiner et al. (2024b).

The first two weeks comprise a continuous simulation with a restart after the first week to set different flux and boundary condition scaling factors for the second week (see Fig. 4a). ICON-ART provides a restart capability, enabling simulations to be paused, checkpointed, and continued from the same meteorological and tracer fields. Simulated mole fractions are written out hourly and interpolated to station locations using inverse distance weighting, with vertical interpolation based on the two nearest levels and horizontal interpolation from the five closest ICON cells. Only daytime values (11:00–16:00 UTC) are then averaged and used for assimilation. Each week is optimized twice using observations from two weeks, except for the first week of the first cycle (see Fig. 4a) which is optimized only once using the observations from the first two weeks.

The second cycle starts with another simulation for the first week but now using the optimized scaling factors returned by CTDAS to propagate the optimized boundary conditions and fluxes forward in time to serve as optimized CO₂ initial conditions for the next week. Only one single CO₂ tracer is simulated in this case as seen in the left part of Fig. 4b. The second week of the second cycle is then restarted from these optimized initial conditions and a new ensemble is generated based on the once optimized scaling factors for week 2 from the previous cycle. Consequently, the spread among ensemble members in week 2 is smaller in the second cycle compared to the first (compare panels a and b). The third week restarts from the checkpoint at the end of the second week, generating a new ensemble with a larger spread. Only observations from the third week are assimilated in this cycle, ensuring no observations are assimilated twice. Finally, panel (c) illustrates the third cycle, where the first two weeks are now fully optimized. This process continues for all subsequent cycles, following the same procedure as outlined for the second cycle.

One design choice for the inversion system involved transferring information between different assimilation cycles. One option would be to use the optimized scaling factors from the previous week as prior values for the next week. Alternatively, one could revert to the original prior values, all set to 1. In the first case, the scaling factors optimized for the previous week are assumed to be a good first guess for the present week. In the second case, the weekly scaling factors are assumed to be independent of each other such that each week should restart from the original prior values (λ=1). We followed the same approach as in (Steiner et al., 2024b), which blends these two extremes: We computed the new priors as a weighted mean with a weight of $\mathrm{1}/\mathrm{3}$ for the original prior and $\mathrm{2}/\mathrm{3}$ for the posterior from the previous cycle. This approach propagates information from the previous cycle while ensuring that, in the absence of new observations, the fluxes gradually relax back to the prior within a few weeks.

Another aspect that could be transferred between optimization cycles is the state vector uncertainty or ensemble spread. However, we decided to keep the uncertainty for a given cycle unaffected by the optimized uncertainties of previous cycles. Transferring the optimized uncertainties to the next week tends to excessively reduce the ensemble spread. Keeping the value unchanged allows the system to remain stable and adapt quickly to sudden changes in flux intensity or background mole fraction errors.

2.3.2 State vector x

To estimate anthropogenic and biospheric CO₂ fluxes separately, the state vector included separate scaling factors for anthropogenic fluxes, gross photosynthetic uptake (GPP), and total ecosystem respiration (RE). Each factor applies to the fluxes within a group of four neighboring grid cells. This grouping reduces the spatial resolution of the inversion but was necessary to save computation time and memory. In addition, eight scaling factors were included to adjust the background CO₂ mole fractions in eight inflow regions following the method outlined in Steiner et al. (2024b). Since an assimilation cycle includes two time windows, each scaling factor is required twice: once for each of the two weeks. For Zurich, the state vector size s is $\mathrm{2}\times (\mathrm{2926}\times \mathrm{3}+\mathrm{8})=\mathrm{17}\mathrm{572}$. For Paris, which has a larger domain with 7680 regions, the size was $\mathrm{2}\times (\mathrm{7680}\times \mathrm{3}+\mathrm{8})=\mathrm{46}\mathrm{126}$ elements.

The state vector elements are scaling factors $\mathit{x}=({\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\lambda }}_{\mathrm{s}})$, so the initial prior state vector x_b simply consisted of ones, ${\mathit{x}}_{\mathrm{b}}=(\mathrm{1},\mathrm{1},\mathrm{\dots },\mathrm{1})$. As mentioned before, from week 3 onwards (see Fig. 4), the prior values were defined by a weighted mean between the initial prior and the posterior λ_i values from the previous week.

Treating GPP and RE as separate components in the state vector is a new approach to urban inverse modeling. Separating these two components is generally difficult, and it is particularly challenging when only daytime observations are assimilated. Nevertheless, the study by Tolk et al. (2011) suggested that optimizing scaling factors for GPP and RE separately performs better than optimizing scaling factors for NEE only even when daytime observations are used. Their conclusion was based on evaluating six different regional CO₂ flux inversion approaches in synthetic model experiments. Similar conclusions were drawn also by White et al. (2019). By estimating GPP and RE separately, we provide the inversion more flexibility to adjust the biogenic fluxes. Optimizing only NEE would have required applying the same scaling factors to GPP and RE, which would preserve any prior errors in the relative magnitudes of the two components. We will present separate results for GPP and RE, but we will also discuss the degree to which they can be separated by analyzing their correlations in the posterior covariances.

2.3.3 Ensemble Square Root Filter approach

In order to optimize fluxes of CO₂, the following Bayesian cost function J of observation-model residuals and differences between prior and posterior fluxes is minimized (Peters et al., 2005):

$$\begin{array}{ll}{\displaystyle }J\left(\mathit{x}\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{x}-{\mathit{x}}^{\mathrm{b}}{)}^{\mathrm{T}}{\mathbf{P}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}^{\mathrm{b}})\\ \text{(1)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathbf{H}\mathit{x}-{\mathit{y}}^{\mathrm{o}}{)}^{\mathrm{T}}{\mathbf{R}}^{-\mathrm{1}}(\mathbf{H}\mathit{x}-{\mathit{y}}^{\mathrm{o}}),\end{array}$$

where y^o denotes the observed mole fractions, H the observation operator (implemented via the ICON-ART model sampled at station locations), x the state vector (flux and background scaling factors), x^b the “background” (or prior) state vector, Hx the simulated mole fraction, P the prior error covariance matrix, and R the observation error covariance matrix. The optimal posterior solution is

$$\begin{array}{}\text{(2)}& {\mathit{x}}^{\text{a}}={\mathit{x}}^{\text{b}}+{\mathbf{PH}}^{\mathrm{T}}{\left({\mathbf{HPH}}^{\mathrm{T}}+\mathbf{R}\right)}^{-\mathrm{1}}\left({\mathit{y}}^{\text{o}}-\mathbf{H}{\mathit{x}}^{\text{b}}\right).\end{array}$$

CTDAS implements an Ensemble Square Root Filter (EnSRF, Whitaker and Hamill, 2002), or an Ensemble Square Root Smoother if more than one window per cycle is used, to solve this problem. A detailed description of the theory and implementation of EnSRFs for inverse modeling is provided by Thanwerdas et al. (2025).

In the EnSRF approach, the matrix P of dimension [s×s] is approximated using an ensemble of N state vectors. Let X denote the ensemble of state vector anomalies of size [s×N], representing deviations from the ensemble mean (e.g., the prior has as mean simply 1). Then P is estimated as

$$\begin{array}{}\text{(3)}& \mathbf{P}={\mathbf{ZZ}}^{\mathrm{T}}\approx {\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\underset{\mathbf{X}}{\underbrace{\mathbf{ZG}}}\underset{{\mathbf{X}}^{\mathrm{T}}}{\underbrace{{\mathbf{G}}^{\mathrm{T}}{\mathbf{Z}}^{\mathrm{T}}}},\end{array}$$

where Z is the lower Cholesky decomposition, i.e., one possible definition of the matrix square root, and G is a Gaussian random matrix of reduced size [s×N] drawn from a standard normal distribution (and ${\mathbf{GG}}^{\mathrm{T}}/(N-\mathrm{1})\approx \mathbf{I}$). This provides a low-rank approximation of the full covariance matrix, such that the problem becomes tractable for high-dimensional inverse problems, as we can now re-write Eq. (2) as

$$\begin{array}{}\text{(4)}& {\mathit{x}}^{\text{a}}\approx {\mathit{x}}^{\text{b}}+{\displaystyle \frac{{\mathbf{XY}}^{\mathrm{T}}}{N-\mathrm{1}}}{\left({\displaystyle \frac{{\mathbf{YY}}^{\mathrm{T}}}{N-\mathrm{1}}}+\mathbf{R}\right)}^{-\mathrm{1}}\left({\mathit{y}}^{\text{o}}-\mathbf{H}{\mathit{x}}^{\text{b}}\right),\end{array}$$

where Y=HX, denotes the ensemble of N modeled mixing ratio deviations from the ensemble mean sampled at the observation locations and times. All these mixing ratios were simulated as separate CO₂ tracers within the same forward ICON-ART run. For our inversions, we chose N=186, following Steiner et al. (2024b).

2.3.4 Observation error covariance matrix R

Uncertainties in the differences between modeled and observed mole fractions are described by the observation error covariance matrix R. It accounts for measurement, transport, and representation errors. Representation errors occur when a model with limited spatial resolution cannot resolve measurements at a single point. In most cases, this matrix is considered to be diagonal, implying statistical independence of model–observation differences. This assumption may be less valid for dense urban networks where stations lie within a few kilometers. While more complex formulations including off-diagonal terms have been proposed (e.g., Ghosh et al., 2021), short sensitivity experiments with such configurations in our system showed only minor effects on flux estimates. Therefore, due to the high computational cost of full EnSRF inversions, we used a diagonal R matrix (see Eq. 5) and assimilated observations serially.

To estimate the diagonal elements, we first calculated centered (i.e., bias-corrected) weekly average root mean square errors (RMSEs) between modeled and observed mole fractions at each station (Fig. 5) and then applied a centered 5 week moving average. The RMSEs were generally larger in Zurich than in Paris because most of the stations were mid-cost sites measuring above rooftop with high sensitivity to surrounding emissions. The largest RMSEs occurred in the period from November to March and the lowest in late summer. Additionally, elevated RMSEs were observed in the first week of April 2023. In Paris, the highest RMSEs were found for the sites gns, cds and jus. This was likely because cds and jus were located within the city, while gns had a data gap during the late summer months, when model-observation differences are typically smaller.

Figure 5 Bias-corrected weekly RMSEs in parts per million (ppm) between modeled and observed CO₂ mole fractions for each station in (a) Zurich (September 2022–September 2023) and (b) Paris (January 2023–January 2024).Background sites in Zurich are indicated in red on the y axis.

The smoothed weekly uncertainties for each station served as the model-data mismatch (MDM) thus generating the R matrix as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathbf{R}=\left[\begin{array}{ccc}{\text{MDM}}_{\mathrm{11}}^{\mathrm{2}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\cdots }& {\text{MDM}}_{nn}^{\mathrm{2}}\end{array}\right],\end{array}$$

where n is the number of observations assimilated in a given cycle.

2.3.5 Prior error covariance matrix P

The prior error covariance matrix P describes uncertainties related to prior flux and background scaling factors, as well as their correlations. These correlations effectively reduce the degrees of freedom, which is necessary to avoid overfitting. The observation networks often lack the density to independently constrain all state vector elements. The correlation length defines the spatial scale of independently resolvable structures and how uncertainties aggregate across the domain.

The structure of P was standardized for both inversions to a common block-diagonal format,

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathbf{P}=\left[\begin{array}{ccc}{f}_{\mathrm{a}}\cdot e^{-\mathbf{D}\oslash {\mathbf{L}}_{\mathrm{a}}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& f_{\mathrm{b}}\cdot e^{-\mathbf{D}\oslash {\mathbf{L}}_{\mathrm{b}}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{P}}_{\text{bg}}\end{array}\right],\end{array}$$

where e^⋅ is implied to be taken element-wise, and we use definitions

• f_a=1, f_b=0.5: the prior uncertainty for the anthropogenic and biospheric categories, respectively. Factors between 0 and 1 correspond to standard deviations of 0 % and 100 % per flux region,

• $\mathbf{D}\in {\mathbb{R}}^{N_{\mathrm{r}}\times N_{\mathrm{r}}}$: distance matrix with entries d_ij between flux regions i and j,

• ${\mathbf{L}}_{\mathrm{a}},{\mathbf{L}}_{\mathrm{b}}\in {\mathbb{R}}^{N_{\mathrm{r}}\times N_{\mathrm{r}}}$: length scale matrices for the anthropogenic and biospheric categories with entries L_ij,a and L_ij,b,

• ⊘: element-wise division,

• ${\mathbf{P}}_{\text{bg}}\in {\mathbb{R}}^{\mathrm{8}\times \mathrm{8}}$: block corresponding to background parameters. It has a banded structure with diagonal elements set to 0.005², corresponding to 0.5 % variance per background inflow region, and off-diagonal correlations between neighboring regions: ±1 and ±2 positions are weighted by factors of 0.5 and 0.2, respectively. The total uncertainty is approximately 1 %.

For Paris, the spatial correlation between flux scaling factors was fixed to length scales of $L_{ij,\mathrm{a}}=L_{\mathrm{a}}=\mathrm{20}$ km and $L_{ij,\mathrm{b}}=L_{\mathrm{b}}=\mathrm{50}$ km for anthropogenic and biospheric fluxes, respectively. The greater length scale for biospheric fluxes arises from the assumption that vegetation responds to environmental drivers in a correlated manner across neighboring regions in the model domain.

In Zurich, spatially varying correlation lengths were considered. Shorter correlation lengths were applied within the city, as the dense observation network can resolve flux gradients. Longer correlation lengths were utilized outside the city. The spatial dependence of L_ij was defined as a linear function of the average distance between flux regions i and j to the city center d_c,ij:

$$\begin{array}{}\text{(7)}& {\displaystyle }{L}_{ij,\mathrm{a}}{\displaystyle }=\left({\displaystyle \frac{d_{\mathrm{c},ij}}{d_{\text{max}}}}\right)L_{\mathrm{a}}+\left(\mathrm{1}-{\displaystyle \frac{d_{\mathrm{c},ij}}{d_{\text{max}}}}\right)L_{\text{urban,a}},\end{array}$$

(and analogously defined L_ij,b) where L_a and L_urban,a are the correlation length limits between which we would like L_ij,a to vary outside and inside the city, and d_max is the maximum distance between flux regions. Using the mean distance d_c,ij was important to preserve the symmetry of the P matrix. As in Paris, we set L_a=20 km for anthropogenic and L_b=50 km for biospheric fluxes, while we used L_urban,a=5 km for anthropogenic and L_urban,b=8 km for biospheric fluxes.

The prior and posterior flux uncertainties presented in the Sect. 3 were computed from the full covariance matrices to account for spatial covariances, similar to the approach described in Steiner et al. (2024b). For each weekly inversion cycle, we calculated:

$$\begin{array}{}\text{(8)}& U=\sqrt{{\mathit{g}}^{\mathrm{T}}\mathbf{P}\mathit{g}},\end{array}$$

where P is the prior or posterior error covariance matrix and g is a vector of 0 and 1 (or a binary mask) corresponding to the selected region of model grid (e.g., selected grid cells inside Zurich, Paris or Île-de-France region). For seasonal or annual means, weekly uncertainties were aggregated using standard error propagation.

3 Results

3.1 Evaluation of the forward simulations compared to observations

3.1.1 European simulation

The main purpose of the European simulation was to provide initial and boundary conditions for the two nested domains. Comparisons with measurements in Europe thus provide information on the magnitude and temporal dynamics of errors in these boundary conditions.

Overall, the mean daily afternoon CO₂ mole fractions from the European simulation showed very good agreement with measurements from the ICOS network. The list of sites included in this comparison is provided in Fig. S2 in the Supplement. Figure 6 presents monthly mean mole fractions, model biases, RMSEs, and correlation coefficients averaged across all sites. The annual mean bias was about 1.24 ppm and the RMSE about 4.19 ppm.

Figure 6 Monthly statistics of mean modeled versus observed afternoon CO₂ mole fractions at ICOS stations in the European model domain. Panel (a) shows the contributions of the different components (background, anthropogenic and biospheric) to total CO₂ mole fractions. Total and background mole fractions are offset by −410 ppm. The remaining panels present (b) monthly mean biases, (c) RMSEs, and (d) Pearson correlation coefficient (r).

The high monthly correlation coefficients (mean value of 0.73) suggest that the model captured most of the day-to-day variability in the observations. However, the error statistics show significant variations between months. For instance, biases in simulated CO₂ mole fractions exhibit a clear seasonal pattern. The largest biases, reaching up to 3 ppm, occurred in August and September. This overestimation is likely related to the challenge in accurately modeling biospheric CO₂ fluxes, which dominate during the warm season. Photosynthetic uptake and respiration processes introduce uncertainties in the prior fluxes that are difficult to constrain at the European scale.

RMSEs were quite stable over the year with a tendency of higher values during the cold season (October–March). This increase may be associated with uncertainties in residential heating, which is one of the main anthropogenic sources. Another important factor that may lead to larger errors during the cold season is stronger vertical stratification and shallower boundary layers.

The lowest correlation of 0.5 occurred in April 2023. This could be related to transition-season effects, where shifts in both biospheric activity (onset of vegetation activity) and heating demand (end of heating period) introduce additional uncertainty. All correlations were calculated as Pearson correlation coefficient r.

3.1.2 Nested simulations over Zurich and Paris

A similar comparison of the prior simulation against measurements was conducted for the two nested simulations using instead the measurement networks in Zurich and Paris (two stations of the Paris network are also ICOS stations). The results are presented in Fig. 7 (blue bars) in terms of statistics per station rather than per month. To be comparable with the results for the European domain, correlation coefficients were first averaged by month and then averaged over all months. Otherwise, high correlation coefficients would be obtained due to the strong seasonal cycle in CO₂, which the model is usually able to capture well. RMSEs and biases of the prior simulation (blue bars) were generally much larger for stations in Zurich than in Paris (see also Table 3). This is due to the different network types with mostly rooftop mid-cost sensors in Zurich as opposed to high-precision measurements on tall towers in Paris. The sites in Zurich are mostly located inside the city and at a much lower altitude above the surface than in Paris. As a result, they are more sensitive to emissions from the city and therefore also to errors in these emissions. For the two background sites, Beromünster and Lägern-Hochwacht in the Zurich domain, the model showed very similar performance as for the sites in the Paris domain.

Figure 7 Daytime (11:00–16:00 UTC) mean bias, RMSE, and Pearson correlation coefficient (r) values of CO₂ mole fractions before and after inversion between ICON-ART and the observations over Zurich (a–c) and Paris (d–f). Blue bars denote the results for the prior simulation, orange bars for the posterior simulation after optimization.

Table 3 Bias, RMSE, and Pearson correlation coefficient for station-averaged prior and posterior CO₂ mixing ratio statistics in the cities of Zurich and Paris.

The biases and RMSEs in the Paris domain were of similar magnitude to the values of the comparison with ICOS sites in the European simulation. The correlation coefficients were also comparable for both simulations, with Paris showing only a slight improvement in the priors, indicating that higher resolution did not substantially increase correlation in this case.

3.2 Inversion results

The inversion significantly improved simulated CO₂ mole fractions compared to observations (orange bars in Fig. 7). Daily RMSEs averaged across all stations decreased in both cities, dropping from 11.83–7.1 ppm in Zurich, and from 5.33–4.25 ppm in Paris. The bias was reduced to near zero, from 2.62–0.05 ppm in Zurich and from 1.36 to −0.11 ppm in Paris. The reduction in the biases was mostly a result of the adjustment of the background scaling factors rather than the flux scaling factors. Finally, the correlation coefficients increased from 0.73–0.83 in Zurich and from 0.76–0.83 in Paris. Although the correlations were already high in the prior simulation, the inversion further improved the temporal agreement between the model and observations. After inversion, the distribution of errors across sites was more uniform, suggesting an improvement in spatial CO₂ gradients.

Figures 8 and 9 show timeseries of weekly prior and posterior flux estimates averaged across Zurich and Paris, as well as their surrounding agglomerations. For Zurich, the larger region encompasses the entire area of the nested simulation. In Paris, it corresponds to the Île-de-France (black contour in Fig. 2).

Figure 8 Time series of prior and posterior CO₂ fluxes in Zurich in (a) the total domain and (b) the city. The remaining panels show the anthropogenic (c) and biospheric (d) CO₂ fluxes in Zurich. Biospheric fluxes are shown in terms of net ecosystem exchange (NEE). Vertical bars denote the 1σ uncertainty in the domain-averaged fluxes computed from the prior and posterior error covariance matrices.

Figure 9 Time series of prior and posterior CO₂ fluxes in Paris in (a) the Île-de-France region and (b) the city. The remaining panels show the anthropogenic (c) and biospheric (d) CO₂ fluxes in the city of Paris. Biospheric fluxes are shown in terms of net ecosystem exchange (NEE). Vertical bars denote the 1σ uncertainty in the domain-averaged fluxes.

In Zurich, most flux adjustments were concentrated within the city, while domain-wide totals remained close to the prior estimates (Fig. 8). Anthropogenic emissions, the dominant source of CO₂ in the urban area, were reduced in almost all weeks, especially during the heating season. The reductions were particularly large in the last weeks of December 2022 and the first week of January 2023. This was a period of comparatively warm temperatures and coincided with the Christmas and New Year period when many residents leave the city for holidays. Our prior estimates of emissions from residential heating accounted for outdoor temperatures by following a heating-degree-days approach (see Brunner et al. (2025) for details). However, the strong reduction in the posterior estimates during warm winter periods suggests that this approach may not accurately capture the influence of temperature. Another factor, which likely contributed to the lower posterior emissions, was the energy crisis driven by the Russo–Ukrainian War. This crisis led to strongly increased prices for gas and electricity. Due to the shortage in primary energy, the Swiss government formulated the goal of a reduction of gas consumption by 15 %, a goal that was exceeded during the heating period October 2022 to March 2023 (https://www.news.admin.ch/de/nsb?id=94439, last access: 18 December 2025). Our prior emissions did not account for this factor. Despite these plausible factors, the reductions in winter were likely too strong as discussed below in Sect. 4. The prior emissions showed a maximum in winter and a minimum in summer, which is expected given the important contribution of heating to total CO₂ emissions in the city. The posterior fluxes showed a very different behavior that does not seem to be realistic, as the posterior fluxes showed excessively strong reductions in late December and early January.

The inversion significantly increased the net uptake of CO₂ by the biosphere in late May/June 2023, suggesting that prior estimates underestimated photosynthetic activity. In other periods, the posterior NEE fluxes remained closer to the prior. In most weeks, NEE was much smaller than the anthropogenic fluxes in the city of Zurich, but in summer, they were of similar magnitude. The increase in posterior uptake in June was partially offset by higher anthropogenic emissions, indicating partial separability between anthropogenic and biospheric fluxes.

Some studies, such as Lauvaux et al. (2016); Lian et al. (2021), focused only on the dormant season to minimize such interference from biospheric fluxes when estimating anthropogenic emissions. Comparing uncertainty reductions between summer and winter thus helps assess how biospheric fluxes influence the inversion constraints. For Zurich, the inversion resulted in notably larger reductions for anthropogenic and total flux uncertainties in winter (both about 90 %) than in summer (about 45 % and 35 %, respectively). This difference may partly result from the inversion attributing model transport errors to emissions during winter shallow boundary layer conditions, as discussed earlier. In summer, Zurich showed relatively higher uncertainty reduction for respiration fluxes (about 33 %) compared to winter (about 20 %), while the uncertainty of photosynthetic uptake by the plants (GPP) remained largely unchanged in both seasons. The uptake uncertainties remaining high are likely due to the lack of assimilated nighttime observations, which would help separate RE and GPP. Detailed information about summer and winter prior and posterior fluxes is provided in the Table S3 in the Supplement.

In Paris, the inversion produced much smaller adjustments compared to Zurich (Fig. 9). The total and anthropogenic CO₂ fluxes remained close to the prior estimates, not only in the Île-de-France region, but also in the city. Biospheric fluxes showed minimal changes in both regions, indicating that the inversion system found little evidence to revise the prior estimates. The posterior uncertainties were significantly reduced for anthropogenic fluxes but remained almost unchanged for biospheric fluxes, suggesting that the observation network was less sensitive to CO₂ exchange with vegetation. This aligns with the fact that NEE fluxes were nearly ten times smaller than anthropogenic fluxes. This may result from the CORINE land cover dataset, which inadequately resolves vegetation in cities.

Overall, the relatively small updates in Paris might be explained by a better initial agreement between the model and observations, as well as by a lower signal-to-noise ratio, with anthropogenic enhancements being less distinct against background variability at the tall towers in Paris compared to the rooftop measurements in Zurich.

Consistent with this, the inversion barely adjusted biogenic flux uncertainties in Paris (e.g., GPP uncertainty reduction was −0.4 % in summer and −6.0 % in winter). This underlines the importance of integrating high-resolution urban vegetation data to better resolve biogenic fluxes within city domains. In Paris, the reductions in total and anthropogenic flux uncertainties were similar in both seasons (about 75 %–78 %), with slightly lower reductions in summer, indicating that anthropogenic emissions were effectively constrained owing to the limited sensitivity to biogenic fluxes.

In both cities, the inversion optimized anthropogenic and biospheric fluxes as well as background CO₂ mole fractions from eight inflow regions (Fig. 10). The Zurich case showed larger changes in background levels and a wider spread between different wind directions. Scaling factors were adjusted by up to 2 %, which corresponds to about 8 ppm given a background of about 400 ppm. Background mole fractions were consistently scaled down when air was advected from the southern sector (SSE and SSW) but scaled up when advected from the north (NNE and NNW). The eastern sectors (ENE and ESE) showed upward corrections in winter but mostly downward corrections in summer.

Figure 10 Optimized scaling factors for background CO₂ mole fractions in (a) Zurich and (b) Paris. The different colors correspond to the 8 different inflow regions.

In contrast, the results for Paris exhibited smaller adjustments, mostly within 1 %, and more consistent between the different wind directions, indicating a systematic overestimation in the prior. For both cities, the adjustments are expected to correct for biases in the boundary conditions provided by the European simulation, which itself is not corrected for biases through data assimilation. The larger corrections for Zurich might be due to the city being more strongly influenced by the European continent due to its position further east from the Atlantic compared to Paris.

Figure 11 summarizes the annual mean CO₂ fluxes and their relative changes for Zurich and Paris. The total fluxes are substantially higher in Paris, reflecting not only its larger area but also its greater population and correspondingly stronger emissions. Net biospheric fluxes (NEE) are only significant for the larger domains (with net uptake in both the prior and posterior fluxes) but are negligible over the two cities. The relative changes are thus driven almost entirely by updates to anthropogenic emissions. In Paris, both domain-wide and city-level fluxes remain close to the prior, though total and anthropogenic fluxes increase by about 7 % in the city. In contrast, more substantial adjustments are found for Zurich: domain-wide fluxes are reduced by approximately 10 %, primarily due to large reductions within the city, where annual anthropogenic emissions decrease by 27 %.

Figure 11 Annual total, anthropogenic and biospheric CO₂ fluxes in Zurich and Paris (a). Relative changes in total and anthropogenic fluxes (b). Uncertainties were propagated from weekly flux uncertainties as described in the Data and Methods section.

The spatial distribution of annual mean anthropogenic CO₂ flux updates for both areas is shown in Fig. 12. Larger adjustments are evident in Zurich, particularly within and upstream of the city, reflecting the high sensitivity of the observation network to these regions. Emissions were reduced in the eastern and western parts of Zurich, corresponding to residential and vegetated areas, with the largest decrease occurring near the western border, an area dominated by vegetation. In contrast, fluxes in Zurich's center were scaled up. In Paris, updates are generally smaller. The results for the Île-de-France show a slight increase north of the city. This is the region covered by the sites gns and cds, which showed particularly large deviations from the ICON-ART model simulation (see Fig. 5). Vegetated zones to the east and west of Paris show minimal changes in anthropogenic fluxes, as expected. Within the city, a minor redistribution of emissions is visible with reductions in the northern districts and increases in the south.

Figure 12 Maps of annual mean anthropogenic CO₂ scaling factor updates in (a) the Zurich model domain, (b) the Paris model domain, (c) the city of Zurich and (d) the city of Paris.

4 Discussion

Posterior anthropogenic emissions in Zurich were likely underestimated (i.e., scaled down too much) in some periods. This was especially evident in late December 2022 and early January 2023. As shown in Fig. S7 in the Supplement for one of the rooftop sites, the ICON-ART model tended to underestimate wind speeds during periods of low winds. This behavior is consistent with previous studies, which showed that mesoscale transport models tend to underestimate wind speed and mixing under stable, low-wind conditions, particularly in urban environments with enhanced surface roughness (Bréon et al., 2015; Lauvaux et al., 2016; Wu et al., 2018). The selected period, from 15 December–13 January, coincides with the largest downward adjustments in emissions. When modeled wind speeds exceeded 1 m s⁻¹, the model captured the observed day-to-day variability reasonably well. However, during calm episodes, the model often simulated wind speeds below 1 m s⁻¹ for several consecutive days, whereas the observations remained higher. These prolonged calm periods result in excessive CO₂ accumulation in the model. The inversion system compensates for this overestimation by reducing the emissions. This highlights how errors in representing the actual meteorology, particularly during low-wind episodes, can introduce artifacts into the inversion.

The network-wide statistics presented in Fig. 13 underscore the systematic nature of these transport errors. The figure shows how prior CO₂ biases and RMSEs vary with modeled wind speed across all sites in Zurich and Paris. In Paris, where sensors are installed on tall towers well above the surface layer, modeled wind speeds fell below 1 m s⁻¹ in less than 2 % of observations. At the rooftop sites in Zurich, in contrast, wind speeds below 1 m s⁻¹ occurred during nearly 15 % of all measurement times. This wind speed threshold marks a sharp transition in model performance. Biases in simulated versus observed CO₂ mole fractions increase to around 9 ppm, and RMSEs exceed 25 ppm. These statistics reflect a regime in which the model systematically underestimates wind speeds allowing CO₂ to accumulate unrealistically. Rather than being isolated events, calm conditions were frequently observed at most sites. The problems of ICON-ART in capturing these situations negatively impact flux estimates during stagnant weather episodes in winter.

Figure 13 Bias and RMSE of CO₂ mole fractions between simulated and observed values as a function of wind speed in Paris and Zurich.

The inversion results revealed significant differences in background CO₂ mole fraction corrections between Zurich and Paris. The larger background CO₂ mole fraction adjustments in Zurich compared to Paris are likely related to the difference in their geographic position. Paris, located closer to the Atlantic Ocean, is exposed to air masses less affected by European anthropogenic emissions. In contrast, Zurich lies further inland and is more influenced by continental sources, resulting in a more complex and variable background signal. This makes it more challenging to differentiate between the influence of local emissions and background changes. Especially in winter when anthropogenic emissions increase not only in Zurich but also across Europe. This is reflected in the background adjustments: during periods of large reductions in anthropogenic fluxes in Zurich, the background CO₂ was actually increased for the NNW, ENE, and ESE wind sectors. This pattern is explained by all three background tower sites (Beromünster, Laegern–Hochwacht, and Birchwil Turm). All of them were showing elevated levels of CO₂ mole fractions higher than those simulated by the model, indicating a regional-scale increase in emissions that was not fully captured in the prior model setup.

Another limitation of the model involves the VPRM model. In our simulations, VPRM relied on satellite indices derived from MODIS observations at a resolution of 250–500 m, which is insufficient for resolving urban vegetation. This poses a challenge when modeling photosynthetic uptake in heterogeneous urban environments. For instance, Ren et al. (2017) conducted a spatio-temporal analysis using Landsat data and field measurements and demonstrated that mixed-pixel effects in urban settings severely compromise estimates of tree density and structural vegetation attributes, especially in areas with impervious surfaces such as roads. Likewise, (Velasco et al., 2016) emphasize that urban greenery is often poorly captured in models, and argue that no method currently exists to directly evaluate the CO₂ uptake by urban vegetation. Together, these studies suggest that current urban biosphere models, which are relying on moderate-to-high-resolution satellite imagery, may systematically underestimate the contribution of small vegetation patches and street trees to total CO₂ fluxes, and misrepresent the net carbon exchange in densely built-up areas. In Paris, this issue was likely exacerbated by using the coarse CORINE land cover dataset at 100 m resolution.

The impact of biospheric flux uncertainties on the estimation of anthropogenic emissions becomes increasingly important during the summer months. In our VPRM simulations for Zurich, prior net ecosystem exchange (NEE) was mostly positive (net source) in the city center, even in summer, with a mean NEE of +0.29 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. In contrast, surrounding vegetated areas functioned as sinks, as expected for the growing season (Fig. 14). The positive NEE in the city center appears unrealistic, suggesting an underestimation of photosynthetic uptake by urban vegetation, likely due to the coarse resolution of satellite data. Respiration in VPRM only depends on temperature but not on satellite indices. The different impacts of low resolution MODIS observations on photosynthesis and respiration may thus explain the unrealistic net positive NEE values.

Figure 14 Net ecosystem exchange (NEE) fluxes in Zurich in summer (June/July/August) 2023. (a) Prior NEE, (b) posterior NEE, and (c) posterior–prior difference.

After optimization, the inversion reduced respiration and enhanced uptake, particularly in surrounding green areas, resulting in a mean NEE of −2.15 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ averaged across the city. This adjustment brings the posterior fluxes into better agreement with local eddy covariance measurements, which show that daytime photosynthetic uptake during summer can exceed concurrent anthropogenic emissions, leading to net negative fluxes in certain areas.

Finally, our inversions used only afternoon CO₂ observations, which primarily constrain daytime fluxes. Assimilating only observations when the boundary layer is fully developed is conventional practice as it minimizes errors introduced by a mismatch between simulated and real vertical mixing, which are largest when the boundary layer is low. This approach is further justified by the fact that nighttime and early-morning fluxes are lower than daytime although not negligible: in Zurich they account for roughly 30 % of the daytime flux, and in Paris for about 30 %–35 % (see Fig. S10 in the Supplement). At the same time, it complicates separation of RE and GPP, as during the night only one of the biospheric flux components is active. Posterior correlations between the two components (0.61 in Zurich, 0.58 in Paris) indicate partial coupling, reflecting that the inversion provides some flexibility to adjust GPP and RE separately. However, these correlations also show that GPP and RE are not fully independent in the inversion. Part of the observational constraint affects both components in a similar way. In this sense, separate optimization mainly prevents errors in the prior relative magnitudes from being preserved, but it does not imply that the system can always distinguish their individual variations. The strength of this separability varies significantly over the year and depends on the available observational information and atmospheric conditions. These limitations are not captured by the analytical posterior uncertainties and will require performing sensitivity tests.

A better separation between anthropogenic and biospheric fluxes might be achieved through assimilation of additional tracers such as NO_x or CO co-emitted with CO₂ or radiocarbon. Radiocarbon and co-emitted species have been measured as part of the ICOS Cities project at the Hardau tower but only for a few months and only using the Eddy covariance technique, which doesn't require very precise absolute calibration. Even if measurements would have been available for longer periods, the co-assimilation of CO₂ with other species comes with its own challenges, either because the species are chemically reactive (like NO_x), because emission ratios between these species and CO₂ are uncertain, or because the measurements are very challenging (as in the case of radiocarbon).

Together, these results illustrate how inversion outcomes are influenced by transport model biases and simplifications in prior flux estimates, which vary by space, season, and flux type. The observed corrections provide a quantitative view of where and when the inversion deviates from the prior, revealing limitations in current models and offering guidance for improving urban CO₂ flux estimation.

Table 4 summarizes the anthropogenic CO₂ emission estimates from this study. It also includes information on the inventory data used as a prior and literature values for Île-de-France and Paris. The inversion reduced prior estimates for Zurich city from 1388.9±156.8–1012.3±38.8 kt yr⁻¹, although no independent literature estimate exists for comparison. For Île-de-France, the prior and posterior values are close to each other, and about 16 % lower than the earlier estimate in (Staufer et al., 2016) between August 2010 and July 2011. In Paris, the posterior estimate of 3580.0±101.9 is slightly higher than the prior 3375±172 kt yr⁻¹ and in good agreement with the estimate for the year of 2020 of 3650±1830 described in (Nalini et al., 2022). Consistent with Lian et al. (2023), the inversion provides more reliable constraints during the dormant period, when biogenic fluxes are minimal. Even though our inversions did not assimilate morning observations as in Lian et al. (2023), the city-scale totals remain in good agreement with both the prior inventory and previous studies, highlighting the robustness of our anthropogenic emission estimates. Overall, the inversion reduced posterior uncertainties compared to priors, demonstrating the added value of atmospheric constraints in emission estimates. We note, however, that the posterior uncertainties in our work only reflect the analytical uncertainty of a Bayesian inversion, which assumes that all errors can be represented by Gaussian distributions with zero mean, i.e. all prior assumptions are unbiased. They do not include other sources of uncertainty such as transport model errors or representation errors introduced by comparing simulations with a bulk land surface scheme, which represents cities only through enhanced surface roughness but without a vertical canopy structure, with observations above rooftops. As a consequence, our reported posterior uncertainties are too low, but without an extensive analysis of potential additional errors, we do not know by how much. Our uncertainties reported for Paris are lower than in previous studies (e.g., Nalini et al., 2022; Lian et al., 2023), which also only accounted for analytical uncertainties but used different inversion approaches, prior uncertainties, and different spatial resolutions.

Table 4 Prior and posterior anthropogenic CO₂ emissions in Zurich, Île-de-France region and Paris based on the results from this study, inventory data and other independent estimates.

5 Conclusions

In this study, we applied the ICON-ART-CTDAS inversion framework at high spatial resolution to estimate urban CO₂ fluxes in Zurich and Paris over a full annual cycle. Simulations were conducted for a central European domain along with nested domains centered on each city. To better capture the influence of complex terrain on mesoscale flow dynamics, the Zurich domain was simulated at a finer resolution (500 m) compared to Paris (1 km).

Forward simulations were evaluated against CO₂ observations in all three domains. Simulated mole fractions for the European domain showed high (monthly mean) correlations with observations from the ICOS network, indicating that the model captures effectively day-to-day variability. Biases and RMSEs were of the order of a few ppm. Additionally, the model reproduced daily variations in wind speed and temperature well (see Supplement).

In Zurich, the inversion reduced anthropogenic emissions by 27 % relative to the prior inventory, resulting in a posterior annual emission of 1012.3±38.8 kt yr⁻¹, whereas in Paris, emissions increased by 7 %, yielding a posterior emission of 3580.0±101.9 kt yr⁻¹. Here we emphasize that the reported posterior uncertainties only reflect the inversion statistics and do not capture other sources of error. They should therefore be considered as lower bounds of the total uncertainty.

The differences in posterior emissions between Zurich and Paris highlight several important lessons for urban inversions. First, inversion results are strongly city-specific: network layout, terrain complexity, and city size lead to distinct flux corrections, background updates, and model–data mismatches. Second, observational network design critically affects inversion sensitivity. Zurich's dense rooftop network, combined with complex terrain, amplifies transport-related biases, particularly during stagnant winter conditions with low wind speeds and shallow mixing layers, which can inflate simulated concentrations and misattribute accumulated CO₂ to local sources. In contrast, Paris tower network and flatter terrain produce smoother regional constraints. These conditions make the inversion less sensitive to both transport errors and biogenic fluxes. Finally, background adjustments also differ between the cities. Paris, located closer to the Atlantic Ocean, is exposed to air masses less affected by European anthropogenic emissions. In contrast, Zurich lies further inland and is more influenced by continental sources, resulting in a more complex signal and larger corrections.

Our results also underscore the importance of accurately representing prior biospheric fluxes in urban areas, especially in summer when vegetation uptake significantly influences total CO₂ fluxes. In Zurich, prior estimates of net ecosystem exchange (NEE) indicated an unrealistic net carbon source in the city center. The inversion corrected this by reducing respiration and enhancing uptake in vegetated areas, aligning posterior fluxes more closely with local eddy covariance observations. This highlights the need for improved representation and parameterization of urban vegetation in biospheric flux models to better capture seasonal and spatial variability. Further studies should explore higher-resolution urban vegetation datasets or urban-adapted models, such as urban VPRM (Hardiman et al., 2017), to better constrain the biogenic contribution. Unfortunately, as pointed out by Stagakis et al. (2025), in-situ observations of urban vegetation are still too limited to thoroughly evaluate biospheric flux models in cities.

In addition, the posterior correlations between GPP and RE scaling factors (0.61 in Zurich, 0.58 in Paris) suggest that these components can be at least partially disentangled with the current inversion framework. Future work should include nighttime CO₂ observations during certain nights where atmospheric boundary layers are not very low to improve the separability of the components. Furthermore, targeted sensitivity tests with perturbed biospheric priors and with optimizing only NEE versus separately optimizing GPP and RE should be performed to further evaluate the robustness of the biogenic flux estimates.

We did not perform meteorological data assimilation or use a meteorological ensemble, limiting our ability to explicitly quantify transport-related uncertainties in the simulated CO₂ mole fractions. To estimate model–data mismatch (MDM), we used a pragmatic approach based on bias-corrected RMSE values calculated for each station. Periods of large model–observation discrepancies, especially under low wind speeds, resulted in exaggerated emission corrections, as the inversion system attempted to compensate for transport biases through flux adjustments. We addressed this issue by rejecting extreme outliers. Although this improved the stability of the inversion, addressing transport related biases more fundamentally will require improved meteorological input, potentially through the use of transport ensembles (e.g., Steiner et al., 2024a; McNorton et al., 2020) or joint optimization of fluxes and meteorology.

As urban CO₂ flux estimation becomes central to net-zero planning, improving transport and biospheric model fidelity and expanding observational constraints will be essential. Projects like ICOS Cities (Lan et al., 2024) have demonstrated the importance of deploying tall-tower eddy covariance, isotopic, and street-level sensors in urban environments to validate emission inventories and improve model constraints on anthropogenic and biospheric fluxes. The CoCO₂-MOSAIC 1.0 dataset (Urraca et al., 2024) is complementary to these observational efforts and offers high-resolution (0.1°) emission priors, which is in line with the crucial role of detailed bottom-up information in urban inversion frameworks.

Overall, our study highlights the sensitivity of urban CO₂ flux inversions to transport model biases, observational network design, and uncertainties in prior flux estimates. Addressing these challenges requires integrated approaches that combine high-resolution meteorology, improved biospheric flux modeling, and expanded observational networks tailored to urban complexity. Future research should explore joint inversion frameworks that incorporate both meteorology and flux uncertainties and leverage transport ensembles to better constrain emissions. These advances will be critical for reliable urban carbon monitoring and for supporting policy efforts aimed at reducing emissions in complex city environments.