Comparing space-based to reported carbon monoxide emission estimates for Europe's iron and steel plants

We use satellite observations of carbon monoxide (CO) to estimate CO emissions from European integrated iron and steel plants, the continent's highest-emitting CO point sources. We perform analytical inversions to estimate emissions from 21 individual plants using observations from the TROPOspheric Monitoring Instrument (TROPOMI) for 2019. As prior emissions, we use values reported by the facilities to the European Pollutant Release and Transfer Register (E-PRTR). These reported emissions vary in estimation methodology, including both measurements and calculations. With the Weather Research and Forecasting (WRF) model, we perform an ensemble of simulations with different transport settings to best replicate the observed emission plumes for each day and site. Comparing the inversion-based emission estimates to the E-PRTR reports, nine of the plants agree within uncertainties. For the remaining plants, we generally find lower emission rates than reported. Our posterior emission estimates are well constrained by the satellite observations (90 % of the plants have averaging kernel sensitivities above 0.7) except for a few low-emitting or coastal sites. We find agreement between our inversion results and emissions we estimate using the cross-sectional flux (CSF) method for the seven most strongly emitting plants, building further confidence in the inversion estimates. Finally, for four plants with large year-to-year variability in reported emission rates or large differences between the reported emission rate and our posterior estimate, we extend our analysis to 2020. We find no evidence in either the observed carbon monoxide concentrations or our inversion results for strong changes in emission rates. This demonstrates how satellites can be used to identify potential uncertainties in reported emissions.

1 Introduction

Integrated iron and steel plants are the highest-emitting point sources of carbon monoxide (CO) in Europe. CO is of particular interest as it is both an important air pollutant and relevant to the greenhouse gas (GHG) balance of the atmosphere as an indirect GHG (Daniel and Solomon, 1998). It is a precursor of ozone and reacts with the cleaning agent OH, thereby increasing the atmospheric lifetime of methane (Jacob, 1999; Wuebbles and Hayhoe, 2002). As a product of incomplete combustion, the majority of CO in our atmosphere is emitted by anthropogenic sources (like road transport and industry) and fires (Zhong et al., 2017). As these combustion processes also emit CO₂, better knowledge of CO can support our understanding of CO₂ emissions (Park et al., 2021; Wu et al., 2022).

The importance of air pollution, both for health effects and for better understanding of our atmosphere, is reflected in regulations by the European Union requiring the reporting of emissions of both GHGs and a large number of air pollutants, including CO, at the facility level (EUR-Lex, 2006). As these reports are an important factor in policymaking, there is a need for verification of these reported emission rates using additional measurements like satellite data. In this study we will use data from the TROPOspheric Monitoring Instrument (TROPOMI) to estimate the emission rate of the 21 highest-CO-emitting European iron and steel plants.

The iron and steel industry has been marked as an important target for de-carbonization, and there has been a push towards near-zero-carbon-emission production of steel (Skoczkowski et al., 2020; Shahabuddin et al., 2023). However, the largest part of planned capacity (increase) still relies on carbon-intensive production and most of the near-zero production projects are currently in testing stages (Higuera and Van Woensel, 2021; Liu et al., 2022). Combined with the continuously increasing global demand for steel, carbon emissions from steel production have roughly doubled from 2000 to 2020 (Bashmakov et al., 2022). Together, the 21 plants considered emit as much CO as Italy, Europe’s fourth-highest CO emitter (E-PRTR, 2023; Denier van der Gon, H. and CoCO2 WP2 team, 2022).

The large amount of released carbon is a result of using coal as a reductor in steel production (Zang et al., 2023). Integrated iron and steel plants cover all processes from iron dust to the production of cast steel. The iron dust contains a lot of oxides, which have to be separated from the iron. However, as the iron dust is too fine to be processed in the blast furnaces, it is first agglomerated during sintering: using hot air and coke, the particle size increases. Due to the low combustion efficiency of fine particles under these conditions (Mohammad et al., 2023), a lot of CO is produced, which is vented into the air with other byproducts (Ho et al., 2013). Subsequently, the combination of sinter and coke is fed into the blast furnace where the oxygen splits from the iron and combines with the carbon molecules of the coke. The liquidized iron is then collected at the bottom of the furnace. Although the reduction of the iron results in a lot of CO, the gas is caught at the top of the furnace and used as fuel (Rackley, 2017). The fraction of carbon in the liquid iron is too high to make steel. Therefore, in a basic oxygen furnace (BOF), oxygen is led through the iron, which binds with the carbon to lower the carbon content to levels appropriate for steel production. Like the blast furnace, the BOF produces a lot of CO that is captured for use as fuel (Rackley, 2017). Annual emission rates for the combination of all processes in the plants can be reported via continuous stack monitoring, extrapolation of discontinuous measurements, or calculations using emission factors in combination with activity and production numbers (E-PRTR, 2023).

Independently of directly measuring emission rates and/or activity-based calculation, emission rates can also be determined based on the resulting CO enhancements in the atmosphere. Previous work on regional emission quantification and analysis using satellite-based concentration measurements has included the use of MOPITT, SCIAMACHY, and TROPOMI (e.g., Gloudemans et al., 2006; Khlystova et al., 2009; Worden et al., 2010; Girach and Nair, 2014; van der Velde et al., 2021). Additionally, the resolution of TROPOMI, down to 7 × 5.5 km² (across × along track), has been shown to be sufficiently high to study individual cities and CO point sources (Tian et al., 2022; Plant et al., 2022; Leguijt et al., 2023; Goudar et al., 2023; Schneising et al., 2024). The coverage of polar-orbiting satellites like TROPOMI allows for consistent investigation of regions all over the world rather than being confined to places with good reporting infrastructure. Even in Europe, this continuous data availability is important as there are some gaps in the data gathered by the European reporting framework. As an example, Slovakia has not reported emissions beyond the year 2017 following a change in reporting format (E-PRTR, 2023). For locations with a continuous record of emissions, we will demonstrate that the satellite data can be used as an independent verification of the reported emissions.

We use CO observations by TROPOMI for 2019 to perform analytical inversions over the largest 21 European point sources of CO to estimate their emission rates and evaluate consistency with reported emissions. In addition, we perform multi-year analyses for sites with large year-to-year differences in reported emission rates and compare our analytical inversions with other satellite-based emission quantification methods.

2 Data and methods

We use TROPOMI carbon monoxide (CO) data in site-specific analytical inversions to estimate annual CO emissions from the 21 largest iron and steel plants in Europe using their emissions as reported to the European Pollutant Release and Transfer Register (E-PRTR) as prior estimates. We will first describe the TROPOMI data product in Sect. 2.1. Section 2.2 covers the prior emission data, and Sect. 2.3 describes the Weather Research and Forecasting (WRF) forward model. In Sect. 2.4 to 2.6, we describe the inversion framework and uncertainty estimation. Finally, in Sect. 2.7, we describe the cross-sectional flux (CSF) method and the concept of wind rotation, which are supplemental methods to analyze emission rates from satellite data.

2.1 TROPOMI carbon monoxide data product

The TROPOMI instrument is a spectrometer on the ESA Sentinel-5 Precursor satellite which flies in a sun-synchronous orbit with an equatorial overpass at 13:30 local time (Veefkind et al., 2012). Its swath of 2600 km allows for daily global coverage at a resolution down to 7 × 5.5 km² (across × along track) for CO. We use the CO operational product version 2.2.0 (Landgraf et al., 2021), which employs the shortwave infrared CO retrieval (SICOR) algorithm to determine the total column CO concentration based on the absorption of reflected sunlight in the shortwave infrared band (SWIR, 2305–2385 nm) (Borsdorff et al., 2018). The ground-based Total Carbon Column Observing Network (TCCON; Wunch et al., 2011) also measures the total column CO concentrations at specific sites by measuring the spectrum of direct sunlight, allowing for validation of the TROPOMI product. TROPOMI shows good agreement with the unscaled TCCON product, with a mean difference per station of 2.45±3.38 % (Sha et al., 2021).

We only use observations with sensitivity to concentrations at the surface and therefore remove observations with a TROPOMI data quality value (QA value) below 0.7 (Landgraf et al., 2021). The remaining data either are cloud-free or contain only low-altitude clouds. Due to the low surface albedo of water, cloudless observations over waterbodies result in more uncertain estimates of the CO concentration. We therefore remove cloudless pixels (pixels with a QA value equal to 1) over water. For all iron and steel plants, we analyze TROPOMI data for 2019. In addition, we analyze 2020 data for four plants, ArcelorMittal Gent (Belgium), ArcelorMittal Gijón (Spain), ArcelorMittal Ostrava (Czech Republic), and Provozovna Třinec (Czech Republic). Because of the heavy computational burden of our analysis, we limit our 2020 analysis to these plants, which show interesting results for 2019 that warrant further investigation.

2.2 Prior emissions: E-PRTR reporting framework and TNO emission inventory

The European Pollutant Release and Transfer Register (E-PRTR) is the official pollutant reporting framework of the European Union (EU) (E-PRTR, 2023). Industries in EU member states are required to annually report facility-level emissions of air pollutants and greenhouse gases to air, water, and soil (EUR-Lex, 2006). For steel production, the reporting requirement applies to all facilities with a capacity exceeding 2.5 t of steel per hour (EUR-Lex, 2006). We use the reported emission rates for 2019, except for U.S. Steel Košice s.r.o. in Slovakia, for which we use the last reported emission rate (2017) instead.

In addition to a reported emission rate, the E-PRTR database contains information on the methods used to determine each specific emission rate as shown in Fig. 1. All measured and calculated emissions obtained conform to either nationally or internationally approved methods. The label “measured” applies to both continuous and short-term discontinuous measurements of the emission rate. “Calculated” emission rates are determined through combined knowledge of activity data (fuel use, steel output) and emission factors, while “estimated” emission rates are determined using non-standardized methods that are not based on publicly available references (ICF, 2020). While Donawitz GesmbH (Austria) mentions the use of stack monitors, which continuously measure the emission rate of gases, and Salzgitter Flachstahl (Germany) reports the use of bi-annual measurements, the majority of the plants do not provide information on the specific method of measurement or calculation, which is in line with findings by ICF (2020).

Figure 1 Facility-level carbon monoxide emissions as reported to E-PRTR from 2017 (top) to 2021 (bottom). The different methods used to determine these emissions are indicated by the colors, measured (green), calculated (orange), and estimated (blue). The year 2019 (middle), which is used for the analysis in this work is shown more opaquely. E-PRTR does not provide uncertainty estimates.

As input to our forward model, we represent anthropogenic CO emissions surrounding the iron and steel plants with the European TNO greenhouse gas and co-emitted species (GHGco) inventory version 4 developed for the EU Horizon CoCO2 project (Kuenen et al., 2022; Denier van der Gon, H. and CoCO2 WP2 team, 2022). The GHGco inventory focuses specifically on Europe and includes emissions for different source sectors grouped following the Gridded Nomenclature For Reporting (Kuenen et al., 2022). A resolution of 0.05° × 0.1° is achieved by combining (1) emission data reported by member states to the Centre on Emission Inventories and Projections of the European Monitoring and Evaluation Programme (EMEP/CEIP), (2) spatial proxies like population density and road networks, and (3) additional datasets like emissions based on reported shipping activity and remotely sensed agricultural fires. The inventory is supplemented with point sources, like iron and steel plants, power plants, and airports, at their exact location. As we use the iron and steel plant emission rate from E-PRTR, we remove the corresponding point sources from the TNO GHGco inventory to avoid double counting of emissions. Because of the CO-intense processes taking place in the iron and steel plants, the emission rates reported to E-PRTR comprise at least 70 % of the total emissions in a 0.4°×0.4° box centered on the plant.

2.3 Forward model: WRF chemical transport model

We use the Weather Research and Forecasting (WRF) chemical transport model version 4.1 (Powers et al., 2017) to simulate three-dimensional concentration fields around each iron and steel plant for 2019. Table A1 shows a list of the plants we use in the simulations and their locations. All simulations use nested domains centered on the location of the plant with an inner domain (147 × 147 km²) simulated at 3 km resolution and an outer domain (441 × 441 km²) at 9 km resolution. Both ThyssenKrupp Schwelgern, and Hüttenwerke Duisburg (Germany) as well as ArcelorMittal Ostrava and Provozovna Třinec (Czech Republic) lie within the inner domain of the other and are combined into one simulation centered midway between the plants.

The E-PRTR emissions for the iron and steel plants are supplemented with anthropogenic emissions from the TNO GHGco inventory (Sect. 2.2). Both the E-PRTR and the TNO emissions are put on a three-dimensional grid using the sector-specific vertical profiles provided by Bieser et al. (2011). The vertical profiles both account for the injection height and include an effective plume rise parametrization. The temporal profiles applied to the emissions per sector are taken from Guevara et al. (2021). Background concentrations are simulated using the 0.25°×0.25° resolution air pollutant forecast product of the Copernicus Atmosphere Monitoring Service (CAMS) as initial and 6-hourly boundary conditions (Inness et al., 2015). In addition to the iron and steel plant, we simulate each sector in the TNO GHGco inventory in the inner domain separately, as well as the CAMS-based background and enhancements in the inner domain originating from emissions in the four quadrants (NE, SE, SW, NW) of the outer domain.

We model carbon monoxide as an inert gas using the contiguous United States (CONUS) physics suite provided in WRF as our baseline setup. Over our small model domain, chemical processes have a small impact on the (long-lived) CO enhancements simulated, while the effect of chemistry outside our domain is included in the CAMS boundary conditions. As will be discussed in Sect. 2.4, it is important for our simulated and observed plumes to have minimal spatial mismatch. However, at the kilometer scale of TROPOMI observations, exact plumes can be difficult to model. A way to minimize the mismatch is by simulating multiple plumes per day using various model settings (Maasakkers et al., 2022a). Therefore, we perform eight simulations for each location using four different planetary boundary layer (PBL) schemes and corresponding surface layer physics and two different driving meteorological fields. The different planetary boundary layer schemes (Mellor–Yamada–Janjić (MYJ) turbulent kinetic energy (TKE), Yonsei University (YSU), eddy-diffusivity mass flux quasi-normal scale elimination (EMF-QNSE), and Mellor–Yamada–Nakanishi–Niino (MYNN) 2.5 level TKE scheme) allow for differences in the vertical distribution and dispersion speed. As driving meteorological fields, we use the National Centers for Environmental Prediction (NCEP, 2000) and the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis products (ERA5) (Hersbach et al., 2020).

To be able to directly compare the simulation to the TROPOMI observations, all simulation output is sampled at the TROPOMI overpass matching footprints of the TROPOMI pixels. The three-dimensional simulation output is converted to a total column by applying the TROPOMI averaging kernel (Landgraf et al., 2021).

2.4 Inversion framework

We use an analytical inversion to estimate posterior emissions as described in Brasseur and Jacob (2017). An advantage of the inversion framework over mass balance approaches is that it more precisely resolves transport in the emission estimation. This comes at the cost of a higher computational load, but it improves the accuracy of the emission estimate and allows the inversion method to be applied to low-coverage situations which would pose challenges to mass balance methods. In addition, the inclusion of more data allows for estimation of smaller emissions over longer time periods.

To estimate emission rates, we optimize the cost function J(x), which is defined as the sum of two parts:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{J}\left(\mathit{x}\right)& ={\left(\mathit{x}-{\mathit{x}}_A\right)}^T{\mathbf{S}}_A^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_A\right)\\ & +\mathit{\gamma }{\left(\mathit{y}-\mathbf{K}\mathit{x}\right)}^T{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{K}\mathit{x}\right).\end{array}\end{array}$$

The first part defines a penalty on a departure of the state vector (x) from the prior inventory emission rates (x_A), weighted by the prior error covariance matrix (S_A). The elements of x correspond to annual emissions from the iron and steel plant, the domain-wide emissions from individual sectors in the TNO GHGco inventory, and emissions from the four quadrants of the outer domain as well as the CAMS-based background concentrations. The second term in Eq. (1) defines a penalty on the difference between observed (y) and simulated (Kx) concentrations, weighted by the observational error covariance matrix S_O. Here K is the Jacobian matrix of the simulation model with respect to changes in emissions. The regularization parameter (γ) is used to avoid overfitting to the TROPOMI observations; its determination is discussed in Sect. 2.5. The optimal posterior solution $\widehat{\mathit{x}}$ which minimizes the cost function is given by

$$\begin{array}{}\text{(2)}& \widehat{\mathit{x}}={\mathit{x}}_A+\mathbf{G}\left(\mathit{y}-\mathbf{K}{\mathit{x}}_A\right).\end{array}$$

Here, G is the gain matrix defined as

$$\begin{array}{}\text{(3)}& \mathbf{G}=\mathit{\gamma }\widehat{\mathbf{S}}{\mathbf{K}}^T{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}},\end{array}$$

with $\widehat{\mathbf{S}}$ being the posterior error covariance matrix:

$$\begin{array}{}\text{(4)}& \widehat{\mathbf{S}}={\left(\mathit{\gamma }{\mathbf{K}}^T{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_A^{-\mathrm{1}}\right)}^{\mathbf{-}\mathbf{1}}.\end{array}$$

With $\widehat{\mathbf{S}}$ we can calculate the averaging kernel of our inversion, which gives the sensitivity of the posterior estimate to the true state:

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

where I represents the identity matrix.

To construct S_A, we assume a diagonal shape and an uncertainty of 20 % for the TNO GHGco inventory, in accordance with the 2σ range of 38 % given in Super et al. (2020). We choose an uncertainty of 10 % for the CAMS background following, for example, Maasakkers et al. (2022b) and Naus et al. (2023) and a 50 % uncertainty applied to the four elements, adjusting for inflow from the outer domain reflecting the high uncertainties associated with long-range transport. As these state vector elements affect many observations, they tend to be well constrained by the observations, and changing their prior uncertainty has limited effect on the outcome of the optimization. To allow for enough flexibility in the inversion, we use an uncertainty of 30 % for emissions from the E-PRTR inventory. However, we test the effect of higher and lower uncertainties in our uncertainty calculation (Sect. 2.6). To construct the observational error covariance matrix S_O, we take the standard deviation of the difference between the simulated concentrations sampled to the TROPOMI footprints and the observations as in Maasakkers et al. (2022a).

As the term (y−Kx) is evaluated for each observation, small mismatches in the exact location of the plume between TROPOMI and the simulation will result in underestimated emissions. This effect can be countered by aggregating the observation to a coarser resolution, in which the simulation and TROPOMI do agree on the position of the plume (Naus et al., 2023). We therefore aggregate TROPOMI observations on a 0.1° grid in our inverse analysis, treating each observation as an independent measurement.

Although aggregation reduces the effect of spatial mismatches between simulation and observation, it is not fully eliminated. Following Maasakkers et al. (2022a), the effect of spatial mismatches can be further mitigated by creating an ensemble of spatially different simulated plumes (Sect. 2.3). For each overpass of TROPOMI, the simulated plume that best matches the observed data is selected for the inversion, further lowering the model-driven spatial mismatch between observation and simulation. We determine which simulation matches the TROPOMI observation best by performing daily inversions with all eight simulation outputs and selecting the simulation with the lowest optimized observational cost (second term in Eq. 1). The different plumes are simulated with two different driving wind fields and four PBL schemes (Sect. 2.3) and further expanded by also selecting the simulated plumes an hour before and after the TROPOMI overpass as in Pandey et al. (2019). Figure 2 shows the eight spatially different simulated plumes at the overpass time as well as the TROPOMI plume observed on the same day. The 16 simulated plumes corresponding to the hour before and after the TROPOMI overpass are not shown. Figure 2i shows the lowest optimized (posterior) observational cost, and this configuration will therefore be used for the optimization for this day. To further limit the contribution of spatial concentration mismatches, we remove days which have the 20 % highest optimized observational cost normalized by the number of pixels. This removes days on which, after aggregation, none of the simulated plumes spatially matched the TROPOMI observed plume well. To ensure representative sampling and avoid structural removal of days without clear observed plume signals, we use the optimized observational cost instead of the prior observational cost. The daily inversions are only used for the selection of the best simulation on each day. Afterwards, the best daily simulations are combined into a single Jacobian and prior vector, and then this is used in an inversion spanning the full year to determine an annual scaling factor for each emission element of the state vector (x). The emission estimates for the iron and steel plants can therefore be directly compared with the annual emission rates reported to E-PRTR.

Figure 2 (a) Concentration as measured by TROPOMI over ArcelorMittal Gent in Belgium (indicated by the x) on 25 February 2019. Panels (b)–(e) and (j)–(m) show different prior simulations using NCEP and ERA5 meteorological data respectively. The variation between the different NCEP and ERA5 simulations is caused by different planetary boundary layer schemes and surface layer physics (as indicated in the titles). Panels (f)–(i) and (n)–(q) show the corresponding posterior concentrations. Out of these posterior simulations, panel (i) shows the lowest observational cost, making it the best simulation for this particular day.

Figure 2 shows a strong southwest to northeast gradient in all prior simulated concentration fields (Fig. 2b–e and j–m) which is not observed in the TROPOMI data. Such strong gradients are not often observed in the simulated data, but days that do have them will negatively impact the accuracy of the inversion result as the over- and underestimates in the simulated concentrations will be compensated for by lowering or increasing respectively the emission rates over affected areas. To reduce the impact of mismatches between the simulated and observed background, we allow our inversion to optimize the background at a daily rather than yearly frequency to prevent biases from aliasing into the emission estimate (Naus et al., 2023). We further split the background into a mean background (which is uniform for all pixels per day) and a deviation from the mean (which is any remaining spatial pattern present in the CAMS-based simulated background). These two parts of the background are added individually to the state vector, yielding two state vector elements per overpass of TROPOMI and giving additional flexibility to the inversion. Figure 2f–i and n–o show that this flexibility results in a reduced spatial gradient in the posterior simulations, better matching the TROPOMI observation. Being derived from CAMS, both the mean daily background and the deviation from the mean are given a 10 % uncertainty. Like the state vector elements for transport from the outer domain, the background is well constrained by the large number of TROPOMI observations, resulting in limited sensitivity to the imposed prior uncertainty. Typical adjustments to both the mean background concentration and its gradient range from 0 %–3 %. However, the daily background gradient state vector element is reduced by up to 30 % to better match observations.

2.5 Regularization parameter

Because of the large number of TROPOMI observations and the assumption of a diagonal observational error covariance matrix, there is a risk of overfitting to the observations. We therefore apply regularization parameter γ to the observational part of the cost function (Eq. 1). To determine an appropriate value for γ, we use the L-curve criterion as described in Hansen (1999). As we are mostly interested in correctly quantifying emissions from the iron and steel plant, we reduce the background contribution to the cost function by scaling the background by the mean difference between simulation and observation over the full year before determining γ. The resulting L curve can be found in Appendix B, from which we conclude that γ=0.1 is appropriate.

2.6 Uncertainty analysis

To evaluate the uncertainty in our posterior emission estimates, we perform an ensemble, varying the relevant parameters in our inversion framework. We report the full spread of this ensemble of inversions for each plant as uncertainty. The range over which each parameter was varied can be found in Appendix C. Figure 3 shows the resulting spread in emission rates related to each varied parameter. Not optimizing the background daily and different aggregation resolutions results in large spreads, exceeding those resulting from the use of different wind products and choices in data filtering. As an alternative to using the observational cost function for selecting the best-matching simulation (Sect. 2.4), we select the simulation based on the highest posterior scaling. We include this (potentially biased high) approach within the uncertainty ensemble to obtain a conservative uncertainty range. Although the choice of the regularization parameter has a small effect on the emission estimate for most plants, it affects a few individual sites more heavily than any of the other variables because they have relatively few observations, and the lower γ values then keep the estimates close to the prior. The posterior estimate is relatively insensitive to variation in the prior, showing that the emission estimates are strongly determined by the TROPOMI observations.

Figure 3 The uncertainty related to each parameter for all plants. The x axis shows the relative deviation from the base posterior estimates. The boxplots show the collection of all investigated plants, with the x axis showing the resulting emission rates normalized by their base posterior estimate. Each boxplot consists of the estimates for each plant for the entire spread in the variable as classified in Appendix C.

2.7 Additional quantification methods

We compare and supplement our inversion approach with two additional methods: the cross-sectional flux (CSF) method and an approach based on oversampling and wind rotation. Both methods solely rely on the CO concentrations measured by TROPOMI and a wind field without incorporating any prior knowledge on emission rates or using simulations of atmospheric concentrations.

The CSF method (Varon et al., 2018) is a “mass balance” emission quantification method that calculates the particle flux at different distances from a source. First, CO enhancements are integrated over cross-sections perpendicular to the plume. Multiplied by an effective wind speed, the full integral over each cross-section gives an emission rate estimate. By repeating this procedure at different distances from the source, an average emission rate corresponding to the observed plume is calculated. The simplicity of the method allows for fast application to different locations at the cost of larger uncertainty and a higher minimum emission threshold compared to methods that rely on large atmospheric transport models. We perform the CSF method as in Leguijt et al. (2023) using the effective wind calibration and 10 m altitude winds from ERA5 (Hersbach et al., 2020). The effective wind speed is a parametrization of the actual wind speed, which aims to account for the effects of turbulence, variation in the vertical wind profile, and plume rise. As uncertainty we report the full range of the same ensemble members as used in Leguijt et al. (2023).

To investigate whether year-to-year variation in inversion-based emission rate estimates is consistent with trends in observed CO concentrations, we also perform a method based on an oversampled wind rotation as in Clarisse et al. (2019). Because of variation in the wind direction, plumes on different days will point in different directions, and oversampling measured concentrations without taking wind-information into account will result in a diffuse enhancement. Valin et al. (2013) and Pommier et al. (2013) showed that the spatially averaged concentrations retain a plume-like shape if the enhancements are rotated around the source location such that the wind points in the same direction. Using the approach as developed in Maasakkers et al. (2022b), we oversample wind-rotated concentration fields and use these as an indication of emission trends rather than a determination of absolute emission rates.

3 Results and discussion

We first discuss the performance of our inversion in Sect. 3.1, followed by a comparison of the satellite-based emission estimates with reported emission rates in Sect. 3.2. To explain any differences between the two, we have extended our analysis to 2020 for some of the investigated plants, for which the results are shown in Sect. 3.4. Section 3.3 explores consistency with the model-independent cross-sectional flux (CSF) method.

3.1 Inversion performance

In Fig. 4a we show the difference between the prior simulation and observations for 2019 over ArcelorMittal Gent (Belgium), gridded at 0.05°. Throughout the domain, excluding five pixels at the center, the concentrations measured by TROPOMI exceed the simulated concentrations. Figure 4c shows the corresponding difference plot for the posterior simulation after optimization of the state vector. This figure shows a reduced bias (1.89 to 0.01 ppb) and absolute bias (2.00 to 0.95 ppb) and a higher correlation (0.80 to 0.83) between simulation and observation, as expected from the optimization. In addition, no spatial patterns are visible in the resulting difference map. The difference between prior and observation can largely be explained by differences in the simulated background. Figure 4b shows the difference between simulation and observation where only the background has been optimized. Although the largest part of the domain shows better agreement with the observations than the prior simulation (mean bias: −0.11; mean absolute bias: 0.97; correlation: 0.83), the simulated concentrations above the iron and steel plant, indicated with the black square, show significant differences. This shows that the inversion framework is sensitive to emissions from the iron and steel plant specifically. The corresponding values for the other plants are shown in Appendix D.

Figure 4 The difference between the prior and posterior simulation and TROPOMI observations over ArcelorMittal Gent (Belgium) for 2019 aggregated at 0.05°. The black squares, with dimensions of 0.15° × 0.15°, are centered on the location of the plant. Panel (a) shows the prior difference, where TROPOMI observations are higher than the simulated concentrations throughout most of the region. Panel (b) shows the impact of optimizing the background, showing a smaller difference, except around the iron and steel plant. In the posterior difference shown in panel (c), there is no clear pattern visible in the difference between the simulation and TROPOMI.

To test whether there is a temporal sampling bias in our method, we investigate the distribution of TROPOMI observations throughout the year per plant. The number of observations shows limited variation over the four quarters of the year, ranging from 19.2 % of the annual number of valid observations (fourth quarter, ArcelorMittal Bremen) to 29.6 % (third quarter, Scunthorpe IS Works).

3.2 TROPOMI emission estimates

Figure 5 shows the TROPOMI-based posterior emission estimates compared to the prior emission rates from E-PRTR. Of 21 posterior emission estimates, 10 agree with E-PRTR within their uncertainty range. The TROPOMI estimates show a correlation of 0.86 with the E-PRTR dataset and are on average 17 % lower, as reflected in Appendix E.

Figure 5 Comparison between our posterior emission estimates (purple) and the emissions reported to E-PRTR for 2019. The colors represent the method used to report the emissions to the E-PRTR framework; these are the same as in Fig. 1 except for U.S. Steel Košice s.r.o., where the gray bar reflects the difference in year between the TROPOMI estimate and the reported emissions. The error bars for the TROPOMI-based estimates show the full spread of the inversion ensemble. The diagonal elements of the averaging kernel of our inversion for each steel plant are shown as hashed bars; these elements are a measure of the extent to which the final emission estimates are constrained by the satellite observations. The value of the inversion averaging kernel corresponds to the fraction of the bar that is shaded.

The hashed bars show the diagonal elements of the averaging kernel of our inversion as given in Eq. (5), which reflect the ability to constrain the iron and steel plant emission estimates based on the satellite observations (Jacob et al., 2016). These elements range between 0 and 1, with 1 indicating that the posterior is fully determined by the observations and values close to 0 resulting in posterior estimates that are mostly determined by the prior. Most of the plants, 67 % (90 %), have inversion averaging kernel values above 0.8 (0.7). This is a result of the large number of TROPOMI observations, with each plant having TROPOMI measurements covering at least part of the simulated plume on 150 to 250 d. ISD Dunaferr (Hungary) forms an exception to the high inversion averaging kernel values, with a diagonal element equal to 0.18. This is also the plant with the lowest prior and posterior emission estimate, resulting in a low sensitivity of the observations to the (small) emissions. With 33 % of the plants having inversion averaging kernel values below 0.8, we are limited to annual emission rate estimates. The 6-monthly emission estimates for these plants would result in averaging kernel values as low as 0.45 and, consequently, emission estimates that are strongly driven by the prior value.

The inversion-based emission estimates for the German plants (ArcelorMittal Bremen, ArcelorMittal Eisenhüttenstadt, AG der Dillinger, Salzgitter Flachstahl, ThyssenKrupp Schwelgern and Hüttenwerke Duisburg) agree within error bars with the emission rates determined by Schneising et al. (2024) using TROPOMI data in a mass balance approach. However, for these six plants, our emission estimates lie at the lower edge of their uncertainty estimates. ArcelorMittal Gent, ArcelorMittal Gijón, Dunkerque, ArcelorMittal Ostrava, and Ilva Taranto show considerably lower posterior emission estimates even though their reported emissions are based on measurements (as indicated by the green color). The same applies to U.S. Steel Košice s.r.o. (Slovakia), although the emission rate was reported for the year 2017. This indicates that our estimates may be conservative as a perfect estimate requires a sufficient spatial match between the modeled and observed plume. Additionally, for ArcelorMittal Gent, ArcelorMittal Gijón, and Ilva Taranto, reported emissions drop by 42 %, 40 %, and 24 % respectively for 2020 compared to 2019 (Fig. 1). We investigate whether the disagreement between our posterior estimates and E-PRTR is persistent in 2020 in Sect. 3.4. Port Talbot S Works and Hüttenwerke Duisburg show large uncertainty ranges. The low values for Talbot originate from the γ=1 and yearly-background optimization ensemble members. As Port Talbot is coastal, a lower number of observations and discrepancies between land and water pixels might create difficulties for the inversion framework. A regularization parameter equal to 1 has a risk of overfitting the observational data and, specifically for Talbot, results in a close-to-zero emission rate estimate. Similarly, not allowing for daily optimization of the background could result in differences between land and water pixels being wrongly interpreted as the effect of the plant. However, the other coastal plants (ArcelorMittal Dunkerque, ArcelorMittal FOS, ArcelorMittal Gijón, Ilva Taranto, and Tata IJmuiden) do not show larger uncertainties compared to inland locations.

The high posterior estimates in the uncertainty range for Hüttenwerke Duisburg come from the ensemble members with γ=1 and prior uncertainty equal to 50 %. Both would allow the inversion to wrongly attribute emissions from the neighboring ThyssenKrupp Schwelgern plant (±10 km difference) to Hüttenwerke Duisburg. However, the corresponding ensemble members for ThyssenKrupp Schwelgern are only 15 %–20 % lower than the base inversion, meaning the summed emission for the two plants is considerably higher than the base inversion in these ensemble members. The fact that the posterior estimate for Hüttenwerke Duisburg specifically is very uncertain shows that the inversion has only limited ability to distinguish between two spatially close point sources.

3.3 Consistency with cross-sectional flux method

In addition to making use of inversions, TROPOMI CO data can be used to determine emission rates using simpler mass balance methods like the cross-sectional flux (CSF) method. Leguijt et al. (2023) found a 100 Gg yr⁻¹ lower limit for the CSF method to be trustworthy when applied to TROPOMI data. Seven of the studied plants – ArcelorMittal Dunkerque, ArcelorMittal Gent, ArcelorMittal Dubrowie, Port Talbot S Works, ThyssenKrupp Schwelgern, Hüttenwerke Duisburg, and U.S. Steel Košice s.r.o. – have reported emission rates above 100 Gg yr⁻¹ and therefore merit the application of the CSF method. The other plants have emission rates lower than what Leguijt et al. (2023) report as the lower limit for the CSF method to be reliable.

Figure 6 shows a comparison between the two methods for the applicable plants. The annual CSF estimates show the mean of the estimates on individual orbits, where the number of suitable orbits ranges from 120 (Port Talbot S Works) to 220 (ArcelorMittal Gent). All estimates agree within uncertainty. The CSF estimates of four out of seven plants are higher than the inversion estimates, which could support the notion that the inversion estimates are conservative. ThyssenKrupp Schwelgern and Hüttenwerke Duisburg lie only 10 km from one another, meaning the enhancements resulting from emissions at these locations can overlap. The CSF method assumes singular point-like sources and is therefore not fully applicable to this situation. This can explain the disagreement between both methods at ThyssenKrupp Schwelgern. Although the reported emission rates for Port Talbot S Works and U.S. Steel Košice s.r.o. are above the 100 Gg yr⁻¹ emission threshold for the CSF method, the inversion-based estimates of 58 and 55 Gg yr⁻¹ fall considerably below this value. When applying the CSF method to these locations, we find estimates of 59 and 34 Gg yr⁻¹, which lie well below the 100 Gg yr⁻¹ threshold. Therefore, while we have limited confidence in the CSF-retrieved emission rates, they do support that the emission rates could be lower than reported.

Figure 6 Comparison of inversion-based TROPOMI emission estimates with estimates from the mass balance cross-sectional flux (CSF) method. Emissions reported to E-PRTR are also shown, using the same colors as in Fig. 5 to represent the reporting method. As the lower limit for the TROPOMI-based CSF method was estimated at 100 Gg yr⁻¹, only the plants with prior or posterior estimates above this value are compared.

3.4 The 2020 analysis

As discussed in Sect. 3.2, we extend our analysis to 2020 for four of the plants that show considerable differences between the prior and posterior for 2019: ArcelorMittal Gent, ArcelorMittal Gijón, ArcelorMittal Ostrava, and Provozovna Třinec. The latter two show little variation in reported emissions from 2019 to 2020 (Fig. 7). On the other hand, both ArcelorMittal Gent and ArcelorMittal Gijón show a sharp drop in reported emissions (42 % and 40 % respectively) from 2019 to 2020, with the 2019 ArcelorMittal Gent reported emission exceeding the mean of the surrounding 4 years by 61 %. Figure 7 also shows the inversion results for these four plants for 2019 and 2020. ArcelorMittal Ostrava shows little difference between the 2019 and 2020 posterior estimates, as expected from the limited variation in reported emission rates. The inversion-based emission rate estimate for Provozovna Třinec increases from 60 (48–91) to 87 (77–102) Gg yr⁻¹, despite having no variation in reported emissions. Both ArcelorMittal Gent and ArcelorMittal Gijón, which have a very different prior for the 2019 simulations than for the 2020 simulations, show much less variation in the posterior estimate than in the reported emissions. The prior emission estimates for 2020 actually agree better with the posterior 2019 estimates for both plants. The 18 % reduction in the posterior estimate for ArcelorMittal Gent also lies within the uncertainty range of our estimate, showing no clear indication of a reduction in emission from 2019 to 2020, contrary to what is suggested by the large difference in reported emissions for those years. ArcelorMittal Gijón shows a 22 % increase in the posterior emission rate as opposed to a decline, although this increase lies within the uncertainty range of the 2019 estimate. Correlations between simulation and observation are similar between 2020 and 2019 – an average posterior (prior) correlation of 0.79 (0.74) in 2020 compared to 0.77 (0.75) in 2019 – indicating comparable inversion performances.

Figure 7 Extended inversion analysis for four plants using 2020 TROPOMI data. Each plant shows four bars, including the reported emission as well as the inversion emission estimate for the years 2019 and 2020. As in Fig. 5, the colors of the reported emissions represent how they were derived and the hashed bars show the diagonal averaging kernel values of our inversion.

To verify the lack of variation observed in our posterior estimates, we perform annual wind rotations (Sect. 2.7). Figure 8 shows the results for the four plants for 2019 and 2020. We estimate the enhancements related to the plant as the difference between 0.1° × 0.2° (width × length) boxes downwind and upwind of the plant. To quantify uncertainties in year-to-year comparisons, we vary the dimensions of up- and downwind boxes simultaneously by up to 30 % and report the full spread. Due to its coastal location, the wind rotations over ArcelorMittal Gijón do not converge in a clear plume and cannot be used to estimate a variation in enhancement. Over ArcelorMittal Gent, we find an enhancement of 2.1 (2.0–2.3) ppb for 2019 and 2.0 (1.9–2.1) ppb for 2020. This 5 % (4 %–11 %) decrease in enhancement is more in line with the 20 % (3 %–37 %) reduction in our posteriors than with the 42 % reduction in reported emissions. For ArcelorMittal Ostrava, we find a decrease in enhancement of 22 % (10 %–27 %). This decrease may partly be attributed to misalignment of the plume and wind direction, which is more prominent in 2020. Over Provozovna Třinec, we see a 15 % (0 %–19 %) increase in enhancements, which agrees well with the reported lack of variation between 2019 and 2020. For both ArcelorMittal Ostrava and Provozovna Třinec, the percentage changes in wind-rotated enhancements are consistent within error bars with the year-to-year variations in our posterior estimates.

Figure 8 TROPOMI wind-rotated averaged concentrations over ArcelorMittal Gent, ArcelorMittal Gijón, ArcelorMittal Ostrava, and Provozovna Třinec for 2019 and 2020 oversampled to 0.01°. The mean CO concentration across the scene has been subtracted from each image to be able to compare the different years. The first panel shows the boxes used to calculate the upwind (white) and downwind (black) concentrations.

4 Conclusions

We performed analytical inversions with 2019 TROPOMI satellite data to determine annual carbon monoxide emission rates for 21 European integrated iron and steel plants. These plants are the highest-emitting CO point sources in Europe. We compared our top-down emission rate estimates to bottom-up emission rates reported to E-PRTR at the facility level. In doing this, we not only evaluated limitations of the satellite-based approach, but also identified outliers pointing at uncertainties in the reported data. The E-PRTR emission rates are used as prior estimates in our inversions. Per site, the inversion uses one of eight simulations with different meteorology for each day to reach optimal spatial agreement between observation and simulation. We allow further freedom in the inversion by optimizing the CAMS-based background on a daily basis, rather than performing an annual scaling.

We find that the posterior estimates for plants with reported emission rates above 50 Gg yr⁻¹ (the majority of the plants) can be constrained with the satellite observations (inversion averaging kernel values above 0.7 for 90 % of the plants). Our emission estimates show a high correlation of 0.86 to the emission rates reported to E-PRTR, with 10 out of 21 plants agreeing within estimated uncertainties. For the remaining 12 plants, we find lower posterior emission rates than reported to E-PRTR, suggesting our method gives conservative emission estimates. For one of our six coastal sites, and for a location with a nearby (within 10 km) point source, we find large differences in posterior estimates for different inversion setups, highlighting that results for these plants and other plants in similar situations should be treated with care.

For the seven plants reporting emission rates above 100 Gg yr⁻¹, we perform an additional emission quantification using the CSF method which does not rely on prior emission information. We find that the CSF-based emission estimates agree with the inversion-based estimates for isolated plants. For two plants in the United Kingdom and Slovakia, the CSF estimates fall below the 100 Gg reported as the lower limit for this method to produce reliable results. However, these low estimates do provide additional confidence in our inversion-based estimates that are also significantly lower than the reported emissions.

We expand our analysis to 2020 for four plants that show large 2019 discrepancies. The inversion estimates for 2019 and 2020 agree with each other, showing the robustness of the method. For ArcelorMittal Gent (Belgium) and ArcelorMittal Gijón (Spain), the reported emission rates for 2020 are 40 % lower than those reported for 2019, while they agree with the 2019 and 2020 inversion estimates, raising questions about the reported emissions for 2019. Comparison of wind-rotated oversampled TROPOMI data for 2019 and 2020 for ArcelorMittal Gent also shows no indication of a large difference in the emission rate between the years. This example shows how these satellite analyses can be used to identify uncertainties in reported emissions. In general, the good agreement between our results and reported emissions indicates that our framework can be used as a measurement-based approach to estimate CO emissions from large steel plants where site-specific measurements are limited or not available.

Appendix A: Iron and steel plant locations

Table A1 Location of the investigated iron and steel plants.

Figure A1 Location of the investigated iron and steel plants. The coordinates are given in Table A1.

Appendix B: Regularization factor determination

Figure B1 Determination of the regularization parameter using the L-curve criterion from Hansen (1999) as discussed in Sect. 2.4. Different values of the regularization parameter are indicated in the graph with the y axis showing the cost related to deviation from the prior and the x axis showing the observational cost. Both costs have been normalized by the number of state vector elements and the number of observations respectively. For low gamma values, emission estimates do not deviate from the prior, which results in large differences between simulation and observation. For values of gamma that are too high, the inversion overfits the observations, resulting in a strong increase in the prior cost. Based on the bend in the L curve, we chose a value of 0.1 for the regularization parameter.

Appendix C: Uncertainty in our inversion estimates

To estimate the uncertainty and sensitivity of our inversion-based emission estimates, we perform an ensemble of inversions varying different parameters. Table C1 contains the full list of ensemble members, which we describe in detail here:

1. For the regularization parameter γ in Appendix B, 0.1 has been established as a suitable value. In the ensemble we include inversions with γ=0.01 and γ=1.

2. We optimize a daily background optimization in the base inversion; in the ensemble we include inversions that only optimize a single scaling of the background.

3. The ensemble includes inversions that use only NCEP or ERA5 data. We also add ensemble members that only use the simulations sampled at the TROPOMI overpass time, instead of using the hour before and after overpass.

4. Whereas the base inversion uses the optimized observational cost to select the daily simulation, the ensemble includes inversions which use the maximal posterior scaling for simulation selection.

5. As members in the ensemble, we include inversions that use different pixel and orbit filtering than the base inversion. On the pixel level, we include inversions with a minimum TROPOMI data quality value (QA value) of 1.0 compared to 0.7 in the base inversion. This strict filtering removes up to 87 % of data compared to the base inversion. Due to the low number of observations, we use a regularization value of 0.5 (as opposed to 0.1) for this specific ensemble member, which was determined in the same way as described in Appendix B. In addition to filtering based on the QA value, we include inversions with a lower maximum allowed scattering layer height. The base inversion uses a value of 5 km, which corresponds to a QA value of 0.7. We lower this value to 0.5 km, which corresponds to a QA value of 1.0. This inversion differs from inversion using the QA value = 1.0 filtering as the QA value imposes additional bounds on the aerosol optical thickness. By not further constraining the optical thickness, about 50 % of TROPOMI observations with QA values of 0.7 and higher are retained.

6. Orbit filtering in the base inversion is done by removing the top 20 % of overpasses with the highest optimized observational error per pixel. In the ensemble, we include both cases where we remove 40 % of overpasses and cases in which we retain all overpasses.

7. Within the ensemble, the resolution to which we aggregate the simulations and TROPOMI observations is increased from 0.1 to 0.15°. We also include inversions in which no aggregation is applied.

8. Within our ensemble, we change the prior by up to 30 %.

9. We also vary the prior uncertainty from the 30 % uncertainty used in the base inversion. Our ensemble members include uncertainties of 10 %, 20 %, 40 %, and 50 %.

Table C1 Full range over which variables were varied in the uncertainty ensemble.

Appendix D: Inversion performance

Table D1 Comparison between TROPOMI observations and the model for the different locations using the prior estimates (left), only correcting the background (middle), and using the posterior estimates (right). The subcolumns show the mean absolute error, the mean error, and the correlation between simulation and TROPOMI observation. The values in parentheses represent the correlations within 0.25° of the plants to focus on the effect of scaling the plants’ emission rates. All errors are shown in parts per billion. ArcelorMittal Ostrava and Provozovna Třinec share the same simulation, as do Hüttenwerke Duisburg and ThyssenKrupp Schwelgern.

Appendix E: Posterior estimates

Figure E1 shows the same data as Fig. 5 in a scatterplot with the same color scheme. The lines show linear regressions between posterior and reported emissions for the different reporting methods (measured, calculated, and estimated) and for the entire set of plants. The full comparison shows a high correlation of 0.86 and a slope of 0.83. The slope smaller than 1 reflects the fact that the TROPOMI-based emission estimates are lower than those reported by the facilities for most plants. Of the different subsets, the reported emissions based on estimation show the biggest deviation from 1 in their slope although the correlation is high due to the very small number of data points.

Figure E1 Different representation of the data shown in Fig. 5, including the correlation between the datasets. The gray line uses the full set of investigated plants, whereas the colored lines correspond to subsets using different E-PRTR reporting techniques.