Towards improving top–down national CO₂ estimation in Europe: potential from expanding the ICOS atmospheric network in Italy

To determine the optimal network design for Italy, we conducted a series of perfect-transport OSSEs using the Lund University Modular Inversion Algorithm (LUMIA) inverse system developed by Monteil and Scholze (2021). These OSSEs were based on a simulated set of synthetic atmospheric observations, derived from a set of arbitrary, but realistic estimates of CO₂ surface fluxes across Europe, which we considered as the 'true' flux. To evaluate the contribution of each additional station, we perform inversions adding the synthetic observations of each station and a set of alternative CO₂ surface fluxes considered as priors. We then assess the reconstruction of the true flux and the posterior uncertainty reduction as shown in figure 1.

Zoom In Zoom Out Reset image size

LUMIA is a Bayesian atmospheric inversion system that optimizes regional biosphere CO₂ fluxes at a 0.25^∘ resolution while prescribing anthropogenic, oceanic, and biomass burning emissions. The system minimizes a cost function by assimilating synthetic CO₂ observations (section 2.6) over a one-year period in 2018, with a one-month spin-up at both ends of the optimization. To account for contributions not captured by regional fluxes, LUMIA estimates boundary conditions offline using either the two-step inversion method of Rödenbeck et al (2009), as applied by Monteil and Scholze (2021), or interpolated mixing ratios from global inversions such as CAMS (section 2.3). Regional transport is handled with precomputed observation 'footprints' from FLEXPART (Pisso et al 2019), while prior flux information is sourced from multiple datasets (section 2.4).

The system employs the Lanczos minimization algorithm (Lanczos 1950), which requires cost function gradients computed via forward and adjoint transport models. During optimization, Lanczos also derives the leading eigenvalues and eigenvectors of the posterior error covariance matrix, enabling posterior uncertainty reconstruction. The full prior and posterior covariance matrices account for variances and covariances, incorporating spatial and temporal error correlations. In LUMIA, the prior flux covariance follows a decaying exponential function. Further details on the system are provided in Monteil and Scholze (2021), and Gómez-Ortiz et al (2025).

The LUMIA inversion system uses pre-computed observation footprints from the FLEXPART Lagrangian transport model (v10.4) (Pisso et al 2019). For 2018, FLEXPART was run backward in time (14 days) for each monitoring station shown in figure 1, releasing 10 000 particles continuously over each one-hour observation period. Footprints were computed by aggregating particle residence times in surface grid cells (below 100 m.a.g.l.). To estimate background CO₂ mixing ratios, the final particle positions (latitude, longitude, time, and height) were used to interpolate CO₂ fields from the Copernicus Atmosphere Monitoring Service (CAMS) global dataset (v18r3) (https://ads.atmosphere.copernicus.eu/datasets/cams-global-greenhouse-gas-inversion).

FLEXPART simulations were driven by meteorological fields from the ECMWF ERA5 reanalysis (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5) at a 3 hourly temporal resolution and a 0.25^∘ horizontal resolution. The regional domain extended from 33° to 73° N; 15° W to 35° E.

In OSSEs, a reference set of CO₂ fluxes is considered the 'true' state and is used to generate synthetic observations that mimic real atmospheric measurements. These observations are produced by applying a transport model to the true fluxes and adding random noise to simulate measurement errors. This controlled setup allows for a direct evaluation of an inversion system by comparing recovered fluxes against the known true fluxes.

The true fluxes used in this study included biospheric, fossil fuel, fire, and ocean fluxes, sourced from different datasets. Biospheric fluxes (NEE) were taken from LPJ-GUESS (Wu et al 2023), which provides simulations at a 0.5^∘ resolution and hourly timescales. Fossil fuel emissions were based on EDGAR v4.3.2 (Gerbig and Koch 2023), with MACC-TNO temporal variations (Denier van der Gon et al 2011) and COFFEE-based extrapolations (Steinbach et al 2011), both accessed via the ICOS Carbon Portal. Fire emissions were obtained from CAMS-GFAS https://ads.atmosphere.copernicus.eu/datasets/cams-global-fire-emissions-gfas, while ocean fluxes were from the Mercator Ocean biogeochemical analysis https://doi.org/10.48670/moi-00015, both produced for the CoCO2 project https://coco2-project.eu/.

Biospheric prior fluxes were represented by the VPRM model (Mahadevan et al 2008, Gerbig and Koch 2020), chosen for its differences from LPJ-GUESS, allowing for a clearer assessment of inversion adjustments. Prior fossil fuel emissions were taken from the TNO GHGco v4 dataset ($0.1^{\circ} \times 0.05^{\circ}$ resolution), which was prescribed rather than optimized, given its well-characterized inventory-based nature. Similarly, fire and ocean fluxes were treated as fixed inputs, as they have a minimal impact on the net terrestrial CO₂ flux in Europe. Spatio-temporal correlations for the prior uncertainties were defined to decay exponentially with a length of 500 km and 30 days for the biosphere fluxes. The length of the correlation used in this study is similar to what is used in other inversion systems in Europe (Thompson et al 2020, Monteil and Scholze 2021, Munassar et al 2023).

To generate the synthetic observations for our OSSEs experiments, we ran the regional transport model forward in time using the true assumed fluxes described in section 2.4. By utilizing these assumed true CO₂ fluxes, we created a set of perfect ('truth') observations, which then were perturbed following the observational error statistics to generate the 'synthetic observations'. Such noise was added to reduce the impact of assuming a perfect atmospheric transport model.

For the control experiment (base network), the FLEXPART model was run using the same station configuration as described in the ObsPack dataset product of the actual (ICOS and non-ICOS) observations. Specifically, we used the same coordinates (latitude and longitude), surface altitude, and sampling height as described in the metadata of the ObsPack dataset (see supplementary, table S1). It is important to note that within the ObsPack product, observations are available at different sampling heights, so we selected only the highest sampling station.

For the Italian candidate stations, the model was run at the station locations described in supplementary information table S2. We assumed a sampling height of 100 m.a.g.l. for stations with low ground altitudes and 2 m.a.g.l. for mountain stations (e.g. if the station will be located at 1498 m.a.s.l., the sampling height will be 1500 m.a.s.l.). For all stations (both base network and new ones) with ground altitudes lower than 1000 m.a.s.l., we ran the model from 13:00 to 18:00 local time. For stations situated over mountainous terrain at altitudes higher than 1000 m.a.s.l., the model was run from 00:00 (midnight) to 05:00 local time. Assimilation of mountain observations at night is a common approach in atmospheric CO₂ flux inversion, given that it favors the subsidence conditions characterizing free-tropospheric concentrations and avoids the need to re-solve daytime up-slope flows (Peters et al 2007, Lin et al 2017). Model-data mismatch errors are explained and shown in supplementary tables S1 and S2.

Our base atmospheric network (control experiment) relies on the European ObsPack dataset (ICOS RI et al 2023) for the year 2018 (supplementary information, table S1). This dataset includes both ICOS and non-ICOS stations. In 2018, four ICOS atmospheric stations were available in Italy: Lampedusa (LMP), Monte Cimone (CMN), Ispra (IPR), and Plateau Rosa (PRS). Potenza was incorporated in 2023; therefore, we added it to the base network to make the experiment more relevant to current conditions, as shown in figure 1 and supplementary information table S2.

For selecting candidate stations for the extended network, we considered 23 locations from Italy's existing monitoring infrastructure, as well as potential future stations. From the current infrastructure, we selected 15 ICOS ecosystem stations, 5 non-ICOS atmospheric stations, and 3 potential future stations (Chieti (CHI), Monte Venda (VND), and Col Margherita (MRG)) as listed in supplementary information table S2.

Our first OSSE experiment (excluding the control run) involved running the LUMIA inversion system 23 times, defined as experiment S01A. For each run in experiment S01A, we added one station at a time from the candidate list (shown in the supplementary information, table S3) to our base European network. This approach allowed us to evaluate which of the 23 selected sites contributed the most to improving the reconstruction of the true flux values and reducing prior flux uncertainties. Given the lack of stations in the southern part of Italy, our selection criteria mainly rely on selecting the sites that bring more information to that area of the country through the evaluation of different statistical metrics explained in section 2.7.

In total, we conducted five experiments, each building on the results of the previous one. In the second experiment (S01B), we repeated the process used in S01A with the remaining 22 stations. We followed the same approach in the third experiment (S01C), where the inversion was run with the remaining 21 candidate stations, and in the fourth experiment (S01D), with 20 candidate stations. This process can be seen as an incremental optimization scheme, where each station that provides the best true flux reconstruction and the greatest uncertainty reduction is added to the network and removed from the candidate list. The fifth experiment, defined as S01, was conducted by running the inversion with all 23 sites together.

The performance of each hypothetical station was evaluated by assessing how well the LUMIA system reconstructs the 'known' true fluxes, defined here as those generated by the LPJ-GUESS model. To measure this performance, we computed several statistical metrics, including mean bias (MB), root mean square error (RMSE), and Pearson's correlation coefficient. Biases were calculated as the differences between the true fluxes and the control. The same approach was applied to the prior and posterior fluxes (e.g. for each scenario of the extended network). These metrics were computed for weekly resampled flux data across various regions (Northern, Central, and Southern Italy) as shown in figure 2.

Zoom In Zoom Out Reset image size

To assess the improvement in the relative reduction of prior biases, we calculated the percentage bias reduction $(\mathit{B_R})$ as:

$$\begin{equation} \mathit{B_R} = 100 \times \left[ 1 - \frac{\left|\mathrm{MB}_{\mathrm{posterior}} \right|}{\left| \mathrm{MB}_{\mathrm{prior}} \right|} \right]\end{equation} \tag{ 1 }$$

Similarly, the percentage uncertainty reduction $(\mathit{U_R})$ was computed as:

$$\begin{equation} \mathit{U_R} = 100 \times \left[ 1 - \frac{\mathrm{SD}_{\mathrm{posterior}}}{\mathrm{SD}_{\mathrm{prior}}} \right]\end{equation} \tag{ 2 }$$

Where prior and posterior represent the standard deviation (SD) of the prior and posterior flux uncertainties, respectively. The percentage of uncertainty reduction accounts for both variances and covariances. Specifically, both the prior and posterior covariance matrices were computed as the product of the emission uncertainties and the correlation matrix, as explained in section 2.4.

Figure 3(a) shows that the current base ICOS network in Italy provides a good constraint on carbon fluxes throughout the country. We observe that prior biases are reduced from −1.18 (RMSE: 2.36) to −0.21 (RMSE: 0.67) TgC week⁻¹, corresponding to an 82.5% reduction in bias. With the addition of the CHI station to the network (obtained from experiment S01A), the bias in the seasonal cycle for Italy is further reduced to −0.10 (RMSE: 0.59) TgC week⁻¹, representing a 92% reduction-an additional 10% improvement over the current network. In terms of flux uncertainty reduction, the CHI site contributes to an overall 29% decrease in uncertainty across Italy (figure 4(a)). This impact is primarily driven by a combination of the spatial extent of its footprint (supplementary information, figure S3) and correlation effects, which are further discussed in section 3.2. The representativeness of the CHI site (footprint) extends from central-eastern Italy, between the Apennine Mountains and the Adriatic coastline, toward the north, improving constraints beyond the influence of the existing northern sites.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A significant improvement in biases is also observed in the northern part of Italy (figure 3(b)) (considering only the base network). The prior biases and RMSEs in Northern Italy decrease from −0.65 (RMSE: 1.40) to −0.04 (RMSE: 0.18) TgC week⁻¹ in the control experiment, which represents a bias reduction of 94% and 87% in RMSE. The addition of the CHI station to the network does not have much influence in the northern region, as the biases remain almost the same. These results were expected in some way, since the current network is concentrated mostly in the northern part of the country, as shown in figure 5(a). In this region, the IT-MBO site contributes the most to flux uncertainty reduction, as shown in figure 4(b). Footprints in this area extends from Monte Bondone, a mountain in the northeastern Italian Alps, toward the northwest, reaching Monte Cimone (supplementary information, figure S3).

Zoom In Zoom Out Reset image size

In central Italy (figure 3(b)), we see that the current network reduces the bias from −0.25 to 0.04 TgC week⁻¹. Even though the posterior biases are closer to zero and indicate a substantial reduction in bias error, we also see a slight overcorrection that shifted the bias to a small positive value. CHI site, also has a strong influence in the central region, where the uncertainty flux reduction reach 29%.

In Southern Italy (figure 3(c)), the current base network reduces prior biases by only 23%, with biases decreasing from −0.27 (RMSE: 0.59) to −0.21 (RMSE: 0.54) TgC week⁻¹. Adding CHI to the network further decreases prior biases from −0.21 (RMSE: 0.54) to −0.14 (RMSE: 0.48), representing a 48% reduction and an additional 25% improvement over the existing network. In this southern region, the CHI site provides the strongest constraint on fluxes compared with other sites, with uncertainty reduction reaching 24% versus 21.6%–22.8% (figure 4(a)).

Given the above results, our selection criteria focused primarily on improving the southern region of Italy, where only the Potenza (POT) station is located. We further explored the potential of the existing Italian network by performing the inversion again with the remaining 22 stations (while keeping the CHI site fixed in the network; experiment S01B), to then select the combination that provided the best constraint on fluxes in the southern region.

Our analysis of experiment S01B identified Lecce (ECO), located in the heel of the Italian Peninsula (as shown in figure 3(h), as the station providing the strongest additional constraint on southern fluxes, not only because it reduces the prior biases from −0.27 to −0.05 TgC week⁻¹, an 82% reduction relative to control experiment (23%), but also because it achieves the largest RMSE reduction (0.35) and the highest correlation (R = 0.72) among all 22 sites (figures 4(g) and (h)). Adding CHI and ECO together translates to a 34% overall improvement, which is an additional 59% better than the existing network.

Another site, apart from the combination of CHI and ECO, that also provides a significant prior bias reduction is CGR (bias = −0.04 TgC week⁻¹). However, this station has a lower RMSE reduction (0.40 TgC week⁻¹) and correlation (R = 0.66) compared to the CHI-ECO sites (Supplementary table S7). CHI-IT-CP2 also leads to a significant bias reduction (−0.07 TgC week⁻¹), but its impact on bias reduction in the southern region is limited to 74%, approximately 8% less than the reduction achieved when the ECO station is selected in combination with CHI.

Adding a third station to the new network, while keeping CHI and ECO fixed, as shown in table S8 (supplementary information), results in only marginal improvements in reducing prior biases in southern Italy. Most station combinations produce posterior biases between −0.01 and −0.13 TgC week⁻¹. We selected IT-SAS as a potential third site, despite its slightly higher bias (−0.10 TgC week⁻¹) compared to other options, because it provides the highest correlation (R = 0.78) and the lowest RMSE (0.31 TgC week⁻¹) among the available sites. While IT-SAS shows a slightly higher bias, its lower RMSE and stronger correlation make it a preferable choice. A lower RMSE indicates smaller and more consistent errors, while a higher correlation suggests better agreement with flux variations over time. This trade-off is relevant since biases can often be corrected more easily than high RMSE, which reflects greater variability. Selecting a site with more stable and predictable error characteristics ultimately strengthens the network's ability to constrain regional fluxes. Other potential sites in the southern region, such as IT-COL and IT-CP2, offer slightly better biases (−0.01 and 0.03 TgC week⁻¹, respectively) but result in a smaller RMSE reduction (0.34) and lower correlation (R = 0.74) compared to the CHI-ECO-IT-SAS combination.

A fourth combination suggests that adding IT-PCM to a network that already includes CHI, ECO, and IT-SAS (supplementary information table S9) further improves the correlation (R = 0.82), RMSE (0.28 TgC week⁻¹), and bias (−0.02 TgC week⁻¹). Adding IT-SAS and IT-PCM as the third and fourth sites results in a 93% improvement relative to the current base network. This represents an additional 12% improvement beyond the 82% gain already achieved with CHI and ECO alone.

We also evaluated the impact of incorporating all 23 hypothetical stations (as shown in figure 5(f) and supplementary information figure S4) into Italy's atmospheric network and found that the bias reduction (0.74%) is similar in magnitude to that achieved when CHI-ECO sites are added (82%). This suggests that simply increasing the number of sites in the network does not necessarily improve the statistical metrics (bias, RMSE, or R) of the posterior fluxes. These findings were, to some extent, expected, given that the prior error correlation length used in the inversion was set to 500 km. In the discussion section (4), we further demonstrate that even reducing the correlation length to 100 km remains insufficient to fully take advantage of all 23 sites in the network. Since atmospheric CO₂ observations are sparse in both Europe and Italy, a prior error correlation length of 500 km is still reasonable, as observations are often far apart. Increasing the number of sites in the atmospheric network would require further reducing the correlation length to fully capture the influence of each observation on estimated fluxes at nearby locations.

We saw in the previous section that the CHI station (figure 4(a)) was one of the sites that provides the best constraint on fluxes across Italy (uncertainty reduction of 29%), central (uncertainty reduction of 28%), and southern region of Italy (uncertainty reduction 24%), with an exception the northern region, where Monte Bondone (IT-MBO) contributes to a 33% reduction in uncertainty when aggregated over that area.

At grid-cell scale, as shown in figure 5, we see around the CHI station, uncertainties reduction is approximately 30%, and gradually decreases to 28%–26% near the Apennines (Maiella ridge, 2700 m), which extends from north to south along the Italian peninsula. These results are somewhat expected, as the sensitivity of the CHI site decreases in this area (supplementary information, figure S3). This underscores the importance of the uncertainty correlation length within the prior error covariance matrix, which enables more extensive constraints on fluxes beyond the immediate vicinity of the station.

Figures 5(c)–(e) show that adding other stations such as ECO, IT-SAS, and IT-PCM to the Italy-based network increases the uncertainty reduction in southern Italy, although the impact is less significant compared to considering the CHI station alone. When all 23 sites are added, the uncertainties in the southern region of Italy improve only slightly. However, additional stations in the northern region lead to a substantial reduction in uncertainty, reaching up to 40% around the Monte Bondone (IT-MBO; figure 2) station, with values gradually decreasing to 37% in its vicinity.