The synthetic observations used in this study were created within the ESA-funded SMARTCARB project to prepare for the upcoming CO2M mission. The dataset has been extensively described and used in previous works (Brunner et al., 2019; Kuhlmann et al., 2019, 2020a, 2021), and is openly available (Kuhlmann et al., 2020b).

The synthetic NO₂ and CO₂ vertical columns are based on atmospheric transport model simulations obtained with the COSMO-GHG model at 1 km by 1 km resolution. The model domain covers parts of Germany, Poland and Czechia for the year 2015. The synthetic data were further averaged to 2 km by 2 km satellite pixels along the 250 km wide swath of simulated CO2M satellite orbits. In this study, we mainly use a constellation setup with two satellites, but we carry out additional tests with one to six satellites. Furthermore, we apply a Gaussian random noise with standard deviation of 1.5 × 10¹⁵ molec./cm² to the NO₂ simulations and 0.5 ppm to the XCO₂ simulations. The CO2M mission requirements indicate that the CO₂ precision shall be better than 0.7 ppm for vegetation scenario at solar zenith angle of 50 degrees and the NO₂ precision better than 1.5 × 10¹⁵ molec./cm² (Meijer et al., 2020). Figure 1 shows an example of XCO₂ observations on a 250 km wide CO2M orbit overlaid with the simulated XCO₂ field. The NO _x and CO₂ emissions used as input in the simulations are summarized in Table 1 for six large point sources and the city of Berlin (Figure 1). Table 1 includes also the annual emissions at 11 UTC, which approximately corresponds to the foreseen satellite overpass time of CO2M over Central Europe.

Beirle et al. (2019), Beirle et al. (2021) introduced the divergence method to estimate the NO _x emissions from TROPOMI NO₂ retrievals. Here we provide an overview of the method applied to synthetic CO2M satellite retrievals of both NO₂ and CO₂ vertical columns. A more comprehensive theoretical discussion is given in the supplementary material, including more details on the different assumptions. The divergence method is based on the continuity equation (Jacob, 1999) at the steady state, where the divergence of vector field F (flux) is defined as the difference between emissions E and sinks S: ∇ ⋅ F = E − S . (1) The flux F is defined as F = (F _x , F _y ) = (V ⋅ u, V ⋅ v), where V is the vertical column density observed by satellite, and u and v are the eastward and northward winds, respectively, at the level of the enhanced concentrations. As discussed by Beirle et al. (2019), the NO _x sink can be calculated from the NO₂ columns as S = LV/τ, where τ is the NO _x lifetime generally assumed as 4 hours (as used also in the SMARTCARB simulations) and L is the constant NO _x -to-NO₂ ratio (typically assumed as 1.32 as in Beirle et al., 2011, 2019). In the follow-up paper by Beirle et al. (2021), the sink term is neglected due to the uncertainties in the assumed NO _x lifetime and only the divergence is analyzed. The divergence method can also be applied to CO₂ but since its lifetime is extremely long (in the order of centuries) as compared to NO _x , the sink term can be neglected as well.

For the flux calculation we use the wind information from the European Centre for Medium-Range Weather Forecasts (ECMWF) next-generation reanalysis ERA5 dataset (Hoffmann et al., 2019) given at 0.1° ×0.1° grid size resolution. Following the approach by Fioletov et al. (2015), we use the mean value from the layers at 900, 950 and 1000 hPa. For each satellite pixel, we take the closest point from the wind grid and then temporally interpolate the wind values to the measurement time.

As the divergence operator is linear, the divergence can be calculated from the mean F _x and F _y fields as reported by Beirle et al. (2019, 2021). However, due to missing values in the data, the averaging and divergence operators are not entirely commutative. In our analysis, we found that the divergence fields are less affected by missing data if the divergence operator is calculated before averaging, and thus this option was used throughout the paper. Supplementary Figure S1 in the supplement shows an example of the CO₂ divergence maps calculated before and after the temporal averaging, with the former option showing less noisy patterns.

The partial derivatives, needed for the divergence ∇ ⋅ F = (∂F _x /∂ _x , ∂F _y /∂ _y ), are calculated using second-order central differences. For data points along the edges, the partial derivatives are calculated using single-sided differences. To adapt the original divergence approach (Beirle et al., 2019) to long-lived gases, such as CO₂, we remove the atmospheric background (e.g., as in Hakkarainen et al., 2016) before calculating the divergence as the flux is not linear with the column V due to the changing wind speed. Thus, before the calculation of the divergence, we derive the XCO₂ anomaly ( X C O 2 a ) as X C O 2 a x , y , t = X C O 2 x , y , t − X C O 2 bg t . (2) Here we define the background X C O 2 bg for each orbit as the median over the area of interest. The definition of the background can be tuned case-by-case. The full theoretical derivation of the divergence method, the background removal and other aspects of the approach are discussed in the supplementary material.

In order to calculate source-specific emissions from the enhancements in the averaged divergence/emission fields, the peak fitting approach (as in Beirle et al., 2021) can be applied by fitting the following function, including a Gaussian and a linear term: P x , y = A 2 π σ x σ y exp − x − x 0 2 2 σ x 2 exp − y − y 0 2 2 σ y 2 + m x x − x 0 + m y y − y 0 + b , (3) where A is the estimated source-specific flux. The variables σ _x and σ _y describe the width of the two-dimentional Gaussian function, and x ₀ and y ₀ indicate the location of the source. The terms m _x , m _y and b define the linear background field. Here we estimate all of the variables mentioned above using the Markov chain Monte Carlo (MCMC) toolbox developed by Laine (2008) (available online at https://mjlaine.github.io/mcmcstat/) using an adaptive Metropolis algorithm (Haario et al., 2001). We also estimate the statistical noise of the averaged divergence field using MCMC. The fitted parameters are estimated as the mean values of the posterior distribution and the fitting uncertainties as the standard deviation.

As an alternative to the divergence method, we estimate the NO _x emissions also by fitting the synthetic observations with the exponentially-modified Gaussian (EMG) function (Beirle et al., 2011). Before fitting, we apply the wind rotation technique (Fioletov et al., 2015) by rotating each pixel around the point source according to the wind direction so that all scenes have wind direction from west to east. The resulting rotated mean field is then integrated along the latitudinal dimension to derive the NO₂ line densities. Those are then fitted with a 1D EMG model M (as in Beirle et al., 2011) as follows: M x = Q × e ∗ G x + B . (4) Here Q is the emission factor (E = Q × u), B is the background (both given in molec./m), u is the effective wind speed, and (e*G)(x) is the convolution between the Gaussian function G and the exponential decay e x = exp − x − X / d , x ≥ X d o w n w i n d , 0 , x < X u p w i n d . (5) The variable X is the distance along the wind direction between the line density peak and the point source. The NO _x lifetime can be obtained by dividing the e-folding distance d with the mean wind speed as τ = d/u. We estimate all together five parameters: Q, B, X, d, and the width of the Gaussian function. For the estimation of the parameters, we use the MCMC approach described in Section 2.3. The uncertainties related to this method have been discussed extensively in the literature (e.g., Beirle et al., 2011; Fioletov et al., 2015; Goldberg et al., 2019).

Since the signal-to-noise ratio for NO₂ is higher than that for CO₂, it is generally easier to estimate NO _x emissions than CO₂ emissions from satellite observations. In addition, NO₂ retrievals are less affected by the presence of clouds and more observations can be acquired. Thus, the CO₂ emissions can also be estimated by scaling the NO _x emissions obtained using the divergence method or the EMG fitting with a source-specific NO _x -to-CO₂ ratio. On the other hand, NO _x emissions depend on assumptions on the lifetime and the conversion from NO₂ to NO _x , that are not needed for CO₂.

Here we derive the NO _x -to-CO₂ emission ratio by calculating the NO₂-to-CO₂ ratios r(x) at multiple transects along matching NO₂ and CO₂ plumes (see Section 3.3). We calculate r(x) using linear fit between NO₂ and CO₂ columns at each transect. We then fit an exponential decay function (similarly to Kuhlmann et al., 2021): r x = r 0 ⁡ exp − x u × τ , (6) where r ₀ is the estimated NO _x -to-CO₂ emission ratio, u is the mean wind speed and τ is the lifetime. As noted by Kuhlmann et al. (2021), it is not always feasible to fit the exponential decay to the data if the number of transect is too small. In those cases, we take into account the value of the ratio at the transects near the source. This method is adapted from the approach proposed by Hakkarainen et al. (2021), where the NO _x -to-CO₂ emission ratios were derived from TROPOMI NO₂ and OCO-2 CO₂ observations at Matimba Power Station in South Africa. In that case, due to the narrow swath of OCO-2 observations, one cross-section per plume was used. If the entire plume is visible from the satellite observations, as for several orbits in the synthetic CO2M observations, we can calculate the ratios at multiple transects along the plume.

The computation of the divergence is sensitive to systematic and random errors (see, e.g., Beirle et al., 2021). In this synthetic case study, particular errors are the single-sounding precision errors added to the data (see Section 2.1), and errors in the estimated effective wind fields. The temporal averaging of the divergence computations partially overcomes the influence of random errors, and the averaging increases the signal-to-noise ratio of the divergence map. However, we can additionally “denoise” the total vertical column density fields using computer vision techniques (Koene et al., 2021). In this study we analyse two denoising methods that are described shortly below.

A relatively simple technique for denoising is a mean filter. Here, we spatially convolve a gridded total vertical column density image V(x, y, t) with an N × N kernel K N 2 whose elements have the value 1/N ². For example, a convolution with the kernel K 9 x , y = 1 9 1 1 1 1 1 1 1 1 1 (7) effectively computes a filtered image based on the average of 3 × 3 regions in V(x, y, t). The expectation is that the mean filter suppresses random noise, varying from pixel to pixel. On the other hand, some signal will be lost in the process. In this paper, we refer to the filter of Eq. 7 as “mean filter 9.” We additionally considered the 5 × 5 kernel K ₂₅ as “mean filter 25.”

An alternative denoising technique we applied exploits the colocated CO₂ and NO₂ observations in the synthetic CO₂ images. It is based on a method called block-matching and 3D filtering, or BM3D (Dabov et al., 2007), which performs around the upper bound of possible performance that a denoising technique can achieve (Levin and Nadler, 2011). In short, BM3D works by denoising an image patch-wise in two steps. First, 3D blocks are formed of similar-looking patches, which are collaboratively denoised using a hard thresholding step, and the collaboratively denoised patches are re-aggregated into an initial estimate of the denoised image. The second step takes place the same way, except that the hard thresholding is replaced with a Wiener filter based on the initially denoised image, with respect to a user-specified denoising level σ _BM3D. We refer to Dabov et al. (2007) and Lebrun (2012) for further details. A notable extension is that we set up BM3D to consider the joint information present in CO₂ and NO₂ images, by normalizing the two images to the same dynamic range and then linearly adding 0.75 times the scaled CO₂ image with 0.25 times the scaled NO₂ image, to form a joint image. The locations of the selected similar-looking patches were established in this joint image; while the denoising took place for both this joint image and an image consisting of 0.75 times the CO₂ image minus 0.25 times the scaled NO₂ image. After adding the two denoised images and rescaling, we obtain a new CO₂ image, while the filtering was helped by the higher signal-to-noise ratio of the NO₂ images during the patch selection and denoising steps. In this paper, we refer to “BM3D σ _BM3D,” for example “BM3D 5,” to indicate denoised images using a specified noise level in the Wiener filter.

Figure 2 illustrates the CO₂ divergence calculated from the CO₂ synthetic observations for the year 2015 simulated by the COSMO-GHG model in an optimal case without presence of noise or background field. We consider model simulations at 11:00 UTC assuming clear-sky conditions gridded at 0.05°×0.05° resolution. The largest point sources (Table 1), such as the individual power stations (Boxberg, Jänschwalde, Lippendorf, Schwarze Pumpe, and Turów) and the city of Berlin are visible as enhancements in the divergence map.

In practice, the situation is more complicated than illustrated in Figure 2 as several aspects affect the divergence calculation. For example, the amount of available data can be reduced due to the limited coverage of satellite observations and the persistence of cloudy conditions. The effect of clouds is more restricting for CO₂, as compared to NO _x , since almost completely clear sky conditions are needed for a successful CO₂ retrieval, while partially cloudy conditions (cloud fraction smaller than 30%) are considered suitable for reliable NO₂ retrievals. In addition, the calculation of the CO₂ divergence requires removing the background (about 400 ppm), with the anthropogenic enhancements in the order of 1 ppm. The highest challenge is posed by the instrument noise.

To mimic the analysis of satellite observations we consider different constellation options from the SMARTCARB dataset. The divergence is calculated for the full year 2015 with data filtering for cloud free conditions. Supplementary Figure S2 in the supplement shows the CO₂ divergence calculated for constellations with one up to six satellites. For the remaining part of this paper, we will use the setup with two satellites shown in Supplementary Figure S2B.

Figure 3 illustrates the divergence calculated with different model setups. Figure 3A shows the NO _x divergence calculated from simulations with added noise (standard deviation 1.5 × 10¹⁵ molec./cm²), without any background removal and with cloud fraction smaller than 0.3. Figure 3B shows the CO₂ divergence based on simulations without presence of background or noise, and with cloud fraction limit 0.01. Both NO _x and CO₂ divergence maps have similar spatial features, with enhancements close to the main emission sources, but the CO₂ fields are noisier and less sharp, as expected due to the longer lifetime and the more restrictive cloud fraction limit which reduces the number of available observations (Supplementary Figure S3 in the supplement). Figure 3C shows the CO₂ divergence map after adding artificial noise (standard deviation 0.5 ppm) to the simulations. Figure 3D shows the CO₂ divergence after adding the simulated background and subsequently removing it by calculating the anomaly as the difference from the background (i.e. the median for each orbit over the area covered in Figure 3). The noise addition has much larger effect to the divergence patterns than the background addition (and removal). Figure 3E combines Figures 3C,D, including both the noise and the background (and its removal). Finally, in Figure 3F the XCO₂ data are denoised by using a mean filter with constant 5-by-5 kernel before removing the background and calculating the CO₂ divergence. After denoising the data, the CO₂ divergence patterns in Figure 3F are similar to Figure 3B and Figure 3D. The effect of denoising is further discussed in Section 3.2. To calculate emission estimates for individual point sources we apply the peak fitting approach described in Section 2.3. For CO₂ we use the divergence map shown in Figure 3F, that includes the denoised data (with background removed) that would be available from a two-satellite constellation. We are able to calculate the CO₂ emission values for all the major sources in the area except for Turów Power Station, which does not appear as a point source in Figure 3F. For NO _x we fit the peaks from the emission map given as the sum of the divergence and the sink terms (Supplementary Figure S4).

Figure 4 shows the comparison between the NO _x and CO₂ emissions estimates derived using peak fitting and the annual mean of the emissions at 11 UTC used as input in the model simulations for each source (orange symbols). The emission estimates are also presented in Table 2. The emission values sit generally close to the 1:1 line, with high correlation (correlation coefficients R = 0.94 and R = 0.97 for NO _x and CO₂, respectively) between our estimates and the 11 UTC emissions used as model input. Some differences can be related to the differences in emissions during different seasons which might be not homogeneously sampled.

To study the effect of the instrumental noise on the emission estimates, we sampled different realisations of the random noise for creating the observation error. The effect of the noise is rather small for the NO _x emission estimates, while it is more pronounced for CO₂, as shown in Figure 3. The different emission estimates (after denoising) are presented as grey symbols in Figure 4. The largest variability due to different random noise can be found for the city of Berlin (emission estimates ranging from 36 to 112 kton/day) and for Schwarze Pumpe power station (5–65 kton/day). Supplementary Figure S5 in the supplement shows two examples of the CO₂ divergence map with two different sets of random noise corresponding to very different outcomes in terms of emission estimates. Supplementary Figure S5A corresponds to CO₂ emission estimates of 36 kton/day for Berlin and 65 kton/day for Schwarze Pumpe, while in Supplementary Figure S5B Schwarze Pumpe CO₂ emissions are much smaller (5 kton/day) and the Berlin CO₂ emissions are about 85 kton/day.

The error bars in Figure 4 correspond to the fitting error (standard deviation of the posterior distribution) calculated using the MCMC sampling. This fitting error is only a statistical error estimate that describes how well the 2D Gaussian model fits the divergence/emission fields. Thus, the fitting error is likely to underestimate the true error as it does not include any systematic component.

The NO _x emissions can be also estimated by using the EMG method described in Section 2.4. This method is suitable for relatively isolated sources and it can be challenging to apply in the presence of strong neighboring sources. Dedicated algorithms have been developed to deal with this issue (Fioletov et al., 2017; Verstraeten et al., 2018; Liu et al., 2022). Here we apply the wind rotation technique around each source (Fioletov et al., 2015) and calculate the NO _x line densities. We obtain a successful fit to estimate the emissions only for the city of Berlin and Turów Power Station, which are both relatively distant from the other sources in the area of study. The resulting NO _x emission estimates for Berlin and Turów are 64.6 ± 3.5 ton/day and 38.5 ± 4.3 ton/day, respectively (Figure 4, black symbols). The estimates based on EMG fitting agree within the uncertainties with the emission estimates obtained from peak fitting from the emission maps. We also estimate the NO _x lifetime for Berlin and Turów as 1.9 ± 0.2 h and 1.7 ± 0.3 h, respectively. These values are about 50% lower than the lifetime used in the SMARTCARB simulations (4 h). The underestimation of NO _x lifetime using the wind rotation and EMG approach was also found by de Foy et al. (2014).

As shown in Figure 3, denoising is an essential step when estimating the emissions from the CO₂ divergence fields. In Figure 4, we used a simple denoising based on a mean filter with constant 5-by-5 kernel. To further analyze the effect of denoising, we test also the BM3D method with various denoising parameters (σ _BM3D = 4, 8, 10, 15, 30) as well as the mean filter with 3-by-3 kernel before the calculation of the CO₂ divergence (Figure 5). As expected, increasing the denoising parameter σ _BM3D in the BM3D method leads to increasingly smoother divergence fields. Overall, the denoising with σ _BM3D = 15 (Figure 5D) produces similar patterns than the mean filter 25 (Figure 3F). The mean filter 9 (with 3-by-3 kernel, Figure 5F) remains quite noisy.

Figure 6 shows the comparison between the CO₂ emission estimates from peak fitting and the emissions at 11 UTC used as input in the model simulations with different setups for denoising. Here we use the same realization of the instrumental noise to test the effect of different denoising approaches. The emission estimates for Berlin show quite a large spread, while for the Schwarze Pumpe power station the spread is relatively smaller excluding one outlier (BM3D with σ _BM3D = 30). The correlation between the CO₂ emissions estimates from peak fitting and the assumed emissions is generally high, but the overestimation of the emissions for Berlin causes the correlation coefficients to become lower when the BM3D denoising method is applied. Figure 6 also show the CO₂ emission estimates in the cases where no noise and no CO₂ background are considered as well as the case where no denoising is applied.

As alternative to the divergence approach, we also derive the CO₂ emissions by converting NO _x into CO₂ emissions, using a source-specific NO _x -to-CO₂ emission ratio calculated as described in Section 2.5. For each emission source listed in Table 1 we plot all the available overpasses from the SMARTCARB dataset and calculate the cross-sectional NO₂-to-CO₂ ratios along each observable plume (Hakkarainen et al., 2019, 2021). The ratios are calculated over several (0.1 degrees wide) transects perpendicular to the plume direction (white lines in Figure 7), located at regular distance downwind from the source (0.025 degrees intervals starting at 0.05 degrees from the source). We then fit the NO₂-to-CO₂ ratios using the exponential decay function (Eq. 6) to derive the NO _x -to-CO₂ emission ratio at the source as described in Section 2.5. Figure 7 shows an example of this approach for the Jänschwalde power station. In general, fitting an exponential decay to noisy data is challenging, and the estimates are highly influenced by the values near the emission source. This could potentially be an issue with real NO₂ satellite observations, if the NO emissions are not yet fully oxidised near the source.

Figure 8 includes the mean NO _x -to-CO₂ emission ratios obtained averaging the results from multiple plumes for each power station. In all cases, the true emission value is within the one-sigma error. The main challenge of this approach is the limited number of plumes available to derive the ratios. The estimation process was the most robust for the Jänschwalde power station, which has the largest CO₂ emissions in the area of study and several (13) detectable plumes available for the calculation of the ratios. On the other hand, we were able to identify only four plumes for Turów power station, which showed large variability, especially due to the low CO₂ signal compared to NO₂. Turów has the largest NO _x -to-CO₂ emission ratio, which appears to be captured by our estimates despite the large variability. We attempted a similar approach also for the city of Berlin, but it was difficult do identify clear CO₂ emission plumes from noisy simulations mainly because Berlin has to be considered as an area source, rather than a point source, and our results remained inconclusive.

Table 2 summarizes the NO _x and CO₂ emission values presented in Figure 3 and the mean NO _x -to-CO₂ emission ratios from Figure 8. The last column, includes the CO₂ emissions when the NO _x emissions obtained from peak fitting are converted into CO₂ emissions using the source-specific NO _x -to-CO₂ emission ratios. Overall, the results showed in Table 2 are mostly in agreement (within the uncertainties) with the emission values used as input in the simulations (Table 1). We note that systematic errors are not taken into account in the uncertainties. If the source-specific NO _x -to-CO₂ emission ratios are obtained with the same dataset as the NO _x emissions, possible biases in the NO₂ observations or the assumed constant NO _x -to-NO₂ ratio cancels out when the estimated ratios are applied to convert the NO _x emissions to CO₂ emissions.