Catalog of NO_x emissions from point sources as derived from the divergence of the NO₂ flux for TROPOMI

We present version 1.0 of a global catalog of NO_x emissions from point sources, derived from TROPOspheric Monitoring Instrument (TROPOMI) measurements of tropospheric NO₂ for 2018–2019. The identification of sources and quantification of emissions are based on the divergence (spatial derivative) of the mean horizontal flux, which is highly sensitive for point sources like power plant exhaust stacks.

The catalog lists 451 locations which could be clearly identified as NO_x point sources by a fully automated algorithm, while ambiguous cases as well as area sources such as megacities are skipped. A total of 242 of these point sources could be automatically matched to power plants. Other NO_x point sources listed in the catalog are metal smelters, cement plants, or industrial areas. The four largest localized NO_x emitters are all coal combustion plants in South Africa. About $\mathrm{1}/\mathrm{4}$ of all detected point sources are located in the Indian subcontinent and are mostly associated with power plants.

The catalog is incomplete, mainly due to persisting gaps in the TROPOMI NO₂ product at some coastlines, inaccurate or complex wind fields in coastal and mountainous regions, and high noise in the divergence maps for high background pollution. The derived emissions are generally too low, lacking a factor of about 2 up to 8 for extreme cases. This strong low bias results from combination of different effects, most of all a strong underestimation of near-surface NO₂ in TROPOMI NO₂ columns.

Still, the catalog has high potential for checking and improving emission inventories, as it provides accurate and independent up-to-date information on the location of sources of NO_x and thus also CO₂.

The catalog of NO_x emissions from point sources is freely available at https://doi.org/10.26050/WDCC/Quant_NOx_TROPOMI (Beirle et al., 2020).

1 Introduction

Nitrogen oxides (${\mathrm{NO}}_x=\mathrm{NO}+{\mathrm{NO}}_{\mathrm{2}}$) are key species in air pollution and tropospheric chemistry (Seinfeld and Pandis, 2006). For the prediction of air quality with regional atmospheric chemistry models, accurate and up-to-date NO_x emissions on high spatial resolution are essential (Bouarar et al., 2019). Such data are often difficult to gain for countries with restrictive information policy. In addition, bottom-up emission inventories take several years to be compiled and are thus generally outdated for countries with quickly developing industrial activities.

Spectrally resolved satellite measurements of solar backscattered radiation allow for the quantification of NO₂ and other trace gases absorbing in the UV–vis spectral range by their characteristic spectral absorption structures (Platt and Stutz, 2008; Richter and Wagner, 2011, and references therein). Tropospheric vertical column densities (TVCDs), i.e., concentrations of NO₂ integrated vertically across the troposphere, can be derived by removing the stratospheric contribution and applying the so-called air mass factor (AMF) that depends on the NO₂ profile shape as well as on viewing geometry, surface albedo, aerosols, and particularly on clouds.

NO₂ TVCDs from satellite measurements provide independent information on the spatial distribution and strength of tropospheric NO₂ levels on a global scale since the mid-1990s, allowing for the identification of NO_x sources and quantification of NO_x emissions (e.g., Leue et al., 2001; Martin et al., 2003; Mijling and van der A, 2012; Martin, 2008; Monks and Beirle, 2011, and references therein).

In October 2017, the TROPOspheric Monitoring Instrument (TROPOMI; Veefkind et al., 2012) was launched as a single payload of ESA's Sentinel-5 Precursor satellite mission. TROPOMI provides NO₂ TVCDs on unprecedented high spatial resolution (7.2 × 3.6 km² until 5 August 2019, 5.6 × 3.6 km² thereafter) and with a high signal-to-noise ratio (Geffen et al., 2020). Single TROPOMI overpasses clearly reveal NO₂ plumes downwind from strong NO_x sources like megacities (Lorente et al., 2019) or large power plants (Beirle et al., 2019). In temporal mean NO₂ TVCDs, however, the high spatial resolution is partly lost due to the averaging over plumes with different directions (related to the variability of atmospheric winds).

Beirle et al. (2019) thus proposed to average NO₂ fluxes F=Vu, i.e., TVCDs multiplied with horizontal wind components. Upscaling NO₂ to NO_x and applying the continuity equation for steady state, this directly allows for the quantification of NO_x emissions from the divergence, i.e., the spatial derivative of the mean NO_x flux:

$$\begin{array}{}\text{(1)}& D:=\mathrm{\nabla }\cdot \mathit{F}=E-S,\end{array}$$

with F being the mean NO_x flux, E the NO_x emissions, and S representing NO_x sinks, i.e., the chemical loss of NO_x.

In Beirle et al. (2019), maps of $D=\mathrm{\nabla }\cdot \mathit{F}$ and $E=D+S$ have been derived and NO_x emissions have been localized and quantified exemplarily for Riyadh, South Africa, and Germany. The sink term S was estimated assuming a constant lifetime of τ=4 h, as derived from the downwind decay of NO₂ for Riyadh (Beirle et al., 2011). For spatially extended sources, like megacities such as Riyadh, S contributes significantly to the derived emissions. For point sources, however, such as the large power plants around Riyadh, emissions are dominated (≈90 %) by the divergence term D directly: point sources show up as distinct peaks in the divergence map, which are much sharper than the corresponding peaks in mean TVCD maps, as the NO_x flux increases abruptly at the NO_x sources, resulting in large values of flux derivatives.

Here we extend this study to the global scale, with a particular focus on point sources. Note that due to TROPOMI's spatial resolution of about 5 km, point sources could be individual facilities but also the merged emissions from industrial areas. Point sources are identified and quantified based on peaks above local background in divergence maps D directly rather than emission maps E, which allows for clearer identification of point source peaks, as well as for the classification of ambiguous cases by artifacts in D.

As the calculation of D from the derivative of mean fluxes requires gridding of TROPOMI data on high spatial resolution, the data processing on a global scale is demanding for I/O operations and working memory. Thus, the analysis is only performed around stationary NO_x sources, which are defined based on the magnitude as well as the temporal variability of TROPOMI TVCDs.

From the derived divergence maps, a catalog of NO_x point sources is extracted by a fully automated algorithm.

The paper is organized as follows: in Sect. 2, the input datasets used in this study are specified. The detailed data processing is explained in Sect. 3. Section 4 presents the NO_x point source catalog. In Sect. 5, the limitations and the potential of the catalog are discussed, followed by an outlook and conclusions.

2 Input data

2.1 Tropospheric NO₂ column densities

The NO_x point source catalog is based on NO₂ TVCDs from TROPOMI for the years 2018–2019, using the offline product (with successively increasing algorithm version from 0.11.0 on 1 January 2018 to 1.3.0 on 31 December 2019), as provided by KNMI/ESA via http://copernicus.eu (last access: 21 June 2021). Details of the TROPOMI tropospheric NO₂ product are given in Geffen et al. (2019) and Geffen et al. (2020).

TROPOMI is flying on a sun-synchronous orbit with a local overpass time of about 13:30. The pixel size at nadir was 7.2 km × 3.6 km initially and even improved to 5.6 km × 3.6 km from 6 August 2019 on (Geffen et al., 2020). TROPOMI provides daily global coverage, resulting in quite good statistics already for annual means.

2.2 Meteorological data

Horizontal wind fields u and v as well as air pressure p and air temperature T are taken from reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF). Wind fields are required for the calculation of fluxes. Pressure and temperature are needed for scaling NO₂ up to NO_x, which involves (a) the reaction rate constant of NO+O₃ as function of T and (b) the conversion of climatological O₃ mixing ratios to concentrations depending on T and p (see Sect 3.4).

Until August 2019, ERA-Interim data are used with a truncation at T255, corresponding to ≈0.7^∘ resolution. Since September 2019, ERA-5 data are used with a truncation at T639, corresponding to ≈0.3^∘ resolution (Hoffmann et al., 2019). For both datasets, a preprocessed dataset was created where the 6-hourly model output (00:00, 06:00, 12:00, and 18:00 UTC) was interpolated to a regular horizontal grid with a resolution of 1^∘. For future processing, sampling will be based on finer temporal and spatial grids.

2.3 Ozone climatology

Ozone mixing ratios, used for the scaling of NO₂ to NO_x, were taken from the Earth System Chemistry integrated Modelling (ESCiMo) project (Jöckel et al., 2016), using the RC1SD-base-10a simulation for the years 2000–2010. The monthly mean climatology was calculated from the model fields sampled online along the overpass time of OMI aboard Aura (which is close to the TROPOMI overpass time) using the MESSy SORBIT submodel (Jöckel et al., 2010). As the divergence is sensitive for the added NO_x at the source, the relevant ${\mathrm{NO}}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio is that close to ground. We thus took O₃ concentrations from the lowest model layer.

2.4 Power plant database

The World Resources Institute provides an open-access Global Power Plant Database (GPPD) (Byers et al., 2019). We use this database (v1.2) in order to automatically identify NO_x point sources corresponding to power plants.

The GPPD lists almost 30 000 power plants of all kinds, including solar, nuclear, and hydro power. For our purpose, we created a subset of those power plants using coal, gas, or oil as primary fuel. In addition, power plants with capacities below 100 MW are skipped. The resulting subset of GPPD comprises 4654 power plants of which 2013, 2265, and 376 use coal, gas, and oil as primary fuel, respectively.

3 Data processing

In this section we describe the data processing step by step. Table 1 summarizes the main steps and also lists similarities and differences to the procedure described in Beirle et al. (2019).

Beirle et al. (2019)

Table 1 Processing settings in this study as compared to Beirle et al. (2019).

3.1 NO₂ selection

For this study, we select TROPOMI tropospheric NO₂ column densities $V_{{\mathrm{NO}}_{\mathrm{2}}}$ for the years 2018–2019 with values of the data quality indicator (“qa value”) above 0.75, as recommended in Geffen et al. (2019), and effective cloud fractions (CF) below 0.3. These selection criteria are the same as in Beirle et al. (2019).

In addition, we skip measurements with solar zenith angle (SZA) above 65^∘ for the calculation of fluxes. This strict criterion removes observations for sun being low and implicitly results in a gradual removal of wintertime measurements for midlatitudes, while in Beirle et al. (2019), winter months have been skipped explicitly for Germany. Wintertime measurements are skipped in order to avoid unfavorable viewing conditions, snow-covered scenes, and stronger interference with aged plumes due to longer lifetimes. Moreover, the SZA restriction allows us to simply parameterize the NO₂ photolysis as a function of the SZA (see Sect. 3.4).

Since large parts of the globe are free from stationary NO_x point sources, in particular oceans, deserts, and forests, the processing focuses on potentially stationary sources. For this purpose, a selection mask is defined (Fig. 1), which is based on magnitude as well as the temporal variability of TROPOMI NO₂ TVCDs. The construction of the selection mask is explained in detail in the Supplement.

Figure 1 (a) Mean tropospheric NO₂ column density for 2019. (b) Mask M for the selection of pixels investigated in this study. The construction of M is described in the Supplement. Boxes indicate the regions as defined in Table 2.

3.2 Grid

In order to have maximum sensitivity for point sources, the TROPOMI observations have to be oversampled, requiring a fine grid resolution of less than 3 km. For each TROPOMI orbit, $V_{{\mathrm{NO}}_{\mathrm{2}}}$ is thus gridded to a regular longitude/latitude grid with 0.025^∘ resolution for 61^∘ S to 61^∘ N. Note that there are a few small NO_x sources north of 61^∘, but due to the strict SZA threshold of 65^∘, the flux statistics would be poor for higher latitudes.

Gridding is done per orbit based on linear 2-D interpolation of TROPOMI pixel centers using the griddata function from the Python module SciPy (Virtanen et al., 2020). This approach allows for fast gridding. In addition, there are no discontinuities at the TROPOMI pixel borders, which would lead to extremely high (positive and negative) values of the derivative.

All missing values (qa < 0.75) as well as the outermost pixels on each side of the TROPOMI swath (i.e., the pixels with the highest viewing zenith angles) are set to not a number (NaN). This is necessary in order to restrict the area of interpolated TVCDs to the actual area covered by measurements.

Hereafter we denote the longitude and latitude dimensions as x and y, respectively, in vector indices as well as in the text.

3.3 Meteorological data

The meteorological datasets u, v, p, and T are extracted from ECMWF input data by linear interpolation in three steps:

1. In vertical dimension to an altitude of 300 m above ground for each ECMWF input dataset with 6-hourly resolution. As emissions from point sources are the focus of this study, we consider wind fields representative for the transport of freshly released NO_x emissions from stacks. The choice of altitude of wind fields is further discussed in Sect. 5.2.3. Only pixels with the mask M ≥ 1 are extracted.

2. In time dimension to the orbit time stamp of each TROPOMI orbit, as given in the orbit filename. The actual TROPOMI overpass lags the orbit time stamp by about 33 and 67 min at 60^∘ S and 60^∘ N at solstice ±7 min for summer/winter. Thus the orbit time stamp reflects the wind conditions for recent plume histories.

3. In horizontal dimensions (latitude, longitude) to the 0.025^∘ grid.

3.4 Upscaling of NO₂ to NO_x

In this study we scale the measured NO₂ TVCD to a NO_x TVCD for each TROPOMI pixel. The conversion factor L is calculated according to the photostationary steady state:

$$\begin{array}{}\text{(2)}& L:={\displaystyle \frac{\left[{\mathrm{NO}}_x\right]}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}}=\mathrm{1}+{\displaystyle \frac{\left[\mathrm{NO}\right]}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}}=\mathrm{1}+{\displaystyle \frac{J}{k\left[{\mathrm{O}}_{\mathrm{3}}\right]}}.\end{array}$$

The impact of volatile organic compounds (VOCs) is neglected here as the focus is put on NO_x point sources and thus generally high NO_x concentrations.

The photolysis frequency of NO₂ J is parameterized as a function of the SZA θ by

$$\begin{array}{}\text{(3)}& J=\mathrm{0.0167}\exp (-\mathrm{0.575}/\cos \mathit{\theta }){\mathrm{s}}^{-\mathrm{1}}\end{array}$$

as proposed by Dickerson et al. (1982). This parameterization is “accurate to about 10 % for mostly sunny conditions” for SZA < 65^∘ (Dickerson et al., 1982).

The reaction rate constant k for the reaction of NO with O₃ is parameterized as a function of temperature (in kelvin) by

$$\begin{array}{}\text{(4)}& k=\mathrm{2.07}\times {\mathrm{10}}^{-\mathrm{12}}\times \exp (-\mathrm{1400}/T),\end{array}$$

following the recommendations from Atkinson et al. (2004) and IUPAC (2013).

Near-surface ozone mixing ratios are taken from a climatology based on the ESCiMo model simulation (Sect. 2.3) and converted into concentrations based on T and p from ECMWF.

The derived values for L represent conditions for near-surface pollution. For background NO_x in the upper troposphere, the partitioning would be shifted towards NO. However, any additive background is automatically removed by the calculation of the divergence. Thus, the partitioning derived for near-surface concentrations is appropriate also for correcting the added column caused by a point source.

Figure 2 displays the ratio of temporal means of NO_x and NO₂ TVCDs. Global mean is 1.35 with a SD of 0.08. Values for Riyadh, South Africa, and Germany are 1.22, 1.36, and 1.41, respectively, in agreement with the value of 1.32 ± 0.26 applied in Beirle et al. (2019), which was based on the number given in Seinfeld and Pandis (2006) for polluted conditions around noontime.

Figure 2 Effective ${\mathrm{NO}}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio, i.e., mean tropospheric NO_x (as derived by assuming photostationary state according to Eq. 2 for each TROPOMI pixel) divided by the mean NO₂ column density for 2018–2019. Note that only cloud-free observations with SZA < 65^∘ are considered; thus wintertime measurements at mid- to high latitudes are skipped, and the expected latitudinal dependency of the ${\mathrm{NO}}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio is suppressed. The spatial variability is a consequence of the dependency of the photostationary state on actinic flux, ozone concentration, and temperature (Eq. 2).

Note that the actual [NO_x] $/$ [NO₂] ratio close to a point source might be different in the case of high NO concentrations causing O₃ titration. In this case, however, the divergence method would detect the emitted NO_x downwind from the source as soon as the NO is converted to NO₂ after mixing with ambient air. This results in a spatial smearing of the peak in the divergence map, leading to broader peaks, but the same integral (and thus emissions) for the peak fitting algorithm (Sect. 3.8.2). For the final budget of NO_x emissions, which are determined from the integrated peaks, the final photostationary state is thus still adequate.

3.5 Gridded fluxes and divergence

From gridded NO_x columns and gridded wind fields, the gridded NO_x flux in both the x and y direction is derived for each TROPOMI orbit. Mean fluxes are calculated for the period 2018–2019, where calm wind conditions (w<2 m s⁻¹) are skipped. For grid pixels with less than 25 measurements, fluxes are set to missing values due to poor statistics. Note that we do not explicitly skip winter months for midlatitudes, as in Beirle et al. (2019) for Germany, but they are removed implicitly by the strict SZA threshold of 65^∘.

From the mean zonal and meridional flux maps, the divergence map $D=\mathrm{\nabla }\cdot \mathit{F}$ is calculated, which is the basis for the identification and quantification of point sources below.

3.6 Lifetime and background corrections

In Beirle et al. (2019), emission maps were derived by adding the sink term $S=V/\mathit{\tau }$ to the divergence map. For this, a constant lifetime of τ=4 h was applied, as derived from OMI data for Riyadh (Beirle et al., 2011). In addition, V was corrected by subtracting the regional background, as the lifetime estimate was derived for freshly released, near-surface pollution, while upper-tropospheric background NO₂ generally has a longer lifetime.

As discussed in Beirle et al. (2019), the inclusion of the sink term has significant impact on area sources; it contributes about 50 % of integrated emissions for the Riyadh urban area. For point sources, however, the emission signal is by far dominated by the divergence term, for instance accounting for 87 % of the emissions from power plant “PP9” in Beirle et al. (2019).

Within this study, we do not correct for the sink term S for the following reasons:

• The NO_x lifetime is expected to be different for the diverse conditions in the considered regions, covering a large variability of temperature, humidity, actinic flux, VOC levels, and NO_x levels. Note that the expected strong dependency of mean lifetime on latitude is again suppressed here due to the selection of SZA < 65^∘.

• NO_x lifetimes simulated from global models cannot resolve the nonlinearities caused by point sources. In addition, emission inventories used as input to model runs generally have a time lag; thus emissions are outdated for quickly developing countries. Consequently, modeling the NO_x lifetime for the investigated point sources is challenging and uncertain.

• Identifying point sources in the divergence map D directly is more immediate than in $E=D+S$, as point sources reveal sharper peaks in D than in E (Fig. 2 in Beirle et al., 2019). In addition, the identification of ambiguous candidates as, e.g., caused by inaccurate wind fields is clearer based on D directly (Sect. 3.8.1).

The resulting low bias of point source emissions caused by the missing lifetime correction is discussed in Sect. 5.2.2.

Since no lifetime correction with $S=V/\mathit{\tau }$ is performed, also the background correction for V, which was performed in Beirle et al. (2019) in order to exclude upper-tropospheric NO_x with longer lifetime and different L, is omitted in the current study. Note that the local background of V, as well as any potential offset due to, e.g., stratospheric correction, would affect the lifetime correction but have no impact on D, as any additive term is lost by the calculation of the derivative.

3.7 AMF correction

Tropospheric column densities of NO₂ are derived from the total slant column by subtracting the stratospheric column and applying the so-called air mass factor (AMF). The AMF can be derived as the sum of height-dependent ”box AMFs”, representing the vertical measurement sensitivity, weighted by the relative profile (Wagner et al., 2007). In the operational TROPOMI NO₂ product, averaging kernels (AKs) are provided that are proportional to box AMFs and allow us to correct the AMF for a different vertical profile (Eskes et al., 2003).

Validation studies report on a general low bias of NO₂ TVCD from TROPOMI of a factor of 2 and more for polluted sites (e.g.,  Verhoelst et al., 2021; Judd et al., 2020), caused by a high biased AMF. Part of this bias seems to be related to the a priori profiles which do not resolve the pollution profiles close to sources. In addition, there are indications that the cloud heights used for the TVCD retrieval are biased low (Compernolle et al., 2021) and the albedo maps used are biased high (Griffin et al., 2019), resulting in biased AKs.

In Beirle et al. (2019), an AMF correction was performed for South Africa and Germany by applying the provided AK to a near-surface profile. In this study, we do not apply such an AMF correction, as the effects of biased input albedo and cloud height cannot be corrected a posteriori based on the provided AKs but require a reprocessing of the TROPOMI NO₂ product.

Consequently, the low bias of TROPOMI TVCDs is directly transferred into the NO_x emissions listed in the catalog, as discussed in detail in Sect. 5.2. The low bias is expected to be improved with an updated NO₂ product, which will then be used for deriving an updated version of the point source catalog.

3.8 Iterative peak fitting

We apply a fully automated iterative peak fitting algorithm in order to detect NO_x point sources. The goal is to identify clear point source peaks in the divergence map, where a robust quantification of emissions is possible, while ambiguous cases are skipped. Thus, the resulting catalog of point sources is incomplete; a detailed discussion on various reasons for missing point sources is given in Sect. 5.1. But the remaining point sources listed in the catalog correspond to actual NO_x sources with high confidence.

In each iteration, the following procedure is executed:

1. The grid pixel with highest value of the divergence is considered the point source candidate.

2. Each candidate is classified into different categories and skipped if ambiguous (Sect. 3.8.1).

3. For a promising candidate, a 2-D Gaussian is fitted to the divergence peak, and successful fits are included in the catalog (Sect. 3.8.2).

4. The candidate is removed from the divergence map before searching for the next highest D value (Sect. 3.8.3).

The iteration stops as soon as the maximum value of the divergence is less than 0.2 µg m⁻² s⁻¹, which is <4 % of the initial maximum value of the first iteration. Below this threshold, almost no further point sources could be detected which meet the quality criteria listed below. For future versions, the availability of several years of TROPOMI data is expected to decrease the noise in divergence maps and will probably allow us to decrease this threshold and investigate additional small NO_x point sources.

3.8.1 Pre-classification of candidates

Point source candidates are iteratively defined as the location of maximum divergence in the global map. Before fitting a Gaussian peak, and quantifying emissions, however, artifacts and ambiguous cases have to be excluded.

For this, a pre-classification is done based on the divergence map 30 km around the candidate. As soon as the candidate is classified as one of the following categories, the pre-classification stops. In Fig. 3, examples for each category are shown.

Figure 3 Clippings of the mean divergence (2019–2020) showing examples for the different categories of point source candidates checked during pre-classification. White indicates missing data. The location of power plants from GPPD as well as cities is indicated by green markers. (a) “Gap” candidate near Jebel Ali/Dubai. (b) “Negative” candidate near Tehran. (c) “Area source” around Baghdad. (d) “Area source” candidate covering several power plants (with Vindhyachal and Anpara being the largest) in India.

1. Gap category. If more than 25 % of grid pixels are missing within 8 km around the candidate, it is classified as “gap”. Gaps were found primarily at sandy coastlines and are caused by persistent cloud coverage above threshold, which is probably an artifact of the coarse resolution of the albedo map used for the cloud retrieval. In the later stage of the iteration, gaps also occur in the vicinity of strong point sources which have been already removed from the divergence map. Figure 3a displays an example for a candidate of the gap category around Dubai, where the global maximum of divergence was found, but missing data do not allow for further quantifications.

2. Negative category. If the minimum divergence within 30 km is negative with an absolute value larger than 50 % of the maximum value, the candidate is “negative”.

   Negative values of the divergence generally indicate NO_x sinks. Thus values of D<0 have to be expected downwind from sources. But absolute values should be far lower than the positive values at the place of emissions, as the NO_x removal is taking place over large distances (at a wind speed of 5 m s⁻¹, a lifetime of 4 h corresponds to an e-folding distance of 72 km).

   High absolute values of negative divergence thus cannot be explained by chemical loss of NO_x (except for very short lifetimes) but indicate an inappropriate simplification of the complex 3-D wind fields by a 2-D wind vector on rather coarse spatial resolution. Also changes of NO_x emissions or wind patterns on temporal scales of some hours (whereas Eq. 1 assumes steady state, i.e., neglects the temporal derivative of V) can cause high negative values of D.

   Figure 3b displays an example of high negative divergence around Tehran. Obviously, the divergence method fails here, even though TVCDs show a hotspot of very high NO₂ around Tehran that can even be spotted in the global map (Fig. 1a). The reason for the noisy divergence is the location of Tehran next to the Alborz mountains, where actual wind patterns are not described appropriately by the low-resolution wind fields.

3. Area source category. If the candidate is neither classified as “gap” nor as “negative”, sections of zonal and meridional means are calculated in order to allow for a quick check of the spatial extent of the divergence peak. For both sections, the full width at half maximum (FWHM) W_{first guess} is determined by checking for the first occurrence of a value below half of the maximum in both directions. For point sources, a FWHM of about 12 km has been reported in Beirle et al. (2019). Values of W_{first guess} above 17 km in both x and y thus indicate a NO_x source, which cannot be categorized as single point source, but an area source, which could be a large city or an extended industrial area with several point sources nearby. Figure 3c and d display two examples for candidates classified as “area source”. Figure 3c shows the city of Baghdad. In Fig. 3d, an “area source” consisting of several point sources close to each other is shown, primarily the coal-fired power plants Vindhyachal (5 MW) and Anpara (2 MW) in India.

3.8.2 Gaussian fit

If a candidate passes all pre-classification checks, a 2-D Gaussian on top of a linear background is fitted to the peak in the divergence map:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}G(x,y)=& {\displaystyle \frac{A}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_x{\mathit{\sigma }}_y}}\exp \left({\displaystyle \frac{-(x-\mathrm{\Delta }x{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_x^{\mathrm{2}}}}\right)\exp \left({\displaystyle \frac{-(y-\mathrm{\Delta }y{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)\\ & +m_{x}x+m_{y}y+b,\end{array}\end{array}$$

with the fit parameters A being the peak integral; σ_x and σ_y being the width of the Gaussian; Δx and Δy being the shift of the peak maximum (relative to the first guess candidate location corresponding to the highest D value); and m_x, m_y, and b describing a linear background.

In contrast to Beirle et al. (2019), no rotation of the peak is allowed in order to make the fit fast and stable and in order to be able to interpret the widths σ_x,y as actual distance in latitudinal and longitudinal dimensions.

The parameters are determined by a least-squares fit of G to the divergence map within 22 km around the candidate. As starting values, σ_x,y are set to W_{first guess}/2.355 for both x and y. But since the fit yields a more robust measure of the peak width than the simple FWHM estimate determined during pre-classification, the candidate is again classified as “area source”, if $W_{\text{fit}}:=\mathrm{2.355}\times {\mathit{\sigma }}_{x,y}$ exceeds 17 km for x or y.

Otherwise, the candidate is considered to be a point source, where the fitted parameter A of the Gaussian peak represents the corresponding point source emissions. However, in order to only keep robust emission estimates in the point source catalog, cases with emissions below 0.03 kg s⁻¹ (which has been derived in Beirle et al., 2019 as the detection limit for optimal conditions) or a relative fit error above 30 % are categorized as “uncertain”.

3.8.3 Candidate removal

The candidate has to be removed from the global map before the next iteration step. Removal is implemented by setting the divergence values around the maximum to NaN. Note that in Beirle et al. (2019) the fitted peaks were subtracted instead. However, this would introduce a highly structured residue in the divergence map, which would create several new artificial candidates for the fully automated peak search algorithm. Removing candidates by setting D to NaN avoids such artificial point sources and prevents any later interferences from fit residues from neighboring sources.

Depending on the classification, the following procedure is applied:

• For point sources, an ellipse with 2$\times {\mathit{\sigma }}_{x,y}$ as a semi-major/minor axis is removed, with σ_x,y from the Gaussian fit.

• For categories “gaps” and “uncertain”, all pixels within 22 km around the maximum are removed.

• For category “negative”, a larger area (30 km around the maximum) is removed, as negative artifacts generally occur not at but next to the sources.

• For area sources, also an area of 30 km around the maximum is removed, if classified based on W_{first guess}. If the area source was classified based on W_fit, an ellipse with 2$\times {\mathit{\sigma }}_{x,y}$ is removed as for point sources.

3.9 Identification of point sources associated with power plants

We perform an automated match of the point source catalog with the combustion power plants listed in GPPD. For each point source, we search for GPPD entries within 5 km radius. In the point source catalog, we add

• the integrated capacity of all power plants within 5 km,

• a complete list of the names of all power plants within 5 km, and

• the primary fuel of the power plant with highest capacity within 5 km.

4 Results

4.1 Candidate classification

The iterative peak fitting algorithm yields 7250 candidates, of which 451 are classified as point sources. For 242 of these point sources, a match with GPPD power plants was found. Table 3 lists the classification statistics for the regions defined in Table 2.

Table 2 Definition of regions used for the regional statistics shown in Table 3 and for regional figures shown in the Supplement.

^* Note that the regions are defined such that all considered pixels are covered by a limited number of figures with similar area as far as feasible. Region names are mostly based on continents. For Asia, regions are labeled after the countries dominating the detected point sources, gaining tangibility while condoning some inaccuracies in actual country borders.

Table 3 Number of candidates and their respective classification found for the regions defined in Table 2. The number of point sources in brackets refers to the point sources associated with power plants.

Figure 4 displays regional maps of color-coded D for selected regions. The respective maps for all regions listed in table 2 are provided in the Supplement in PDF format, allowing for loss-free zooming.

Figure 4 Divergence D (color coded) and point source classification (symbols) for (a) the Middle East, (b) India, (c) Europe, and (d) Ukraine/western Russia. Triangles display the point sources listed in the catalog, where matches to GPPD power plants are indicated in magenta. Respective maps for all regions defined in Table 2 are provided in the Supplement (Figs. S2–S14). Also for sake of clarity, non-point sources are only shown for candidates with D_max>1 µg m⁻² s⁻¹ or the integrated divergence within 30 km exceeding 0.1 kg s⁻¹.

The candidate classifications are indicated by symbols, where point sources with/without a power plant match are displayed as triangles in magenta/dark grey, respectively. Non-point sources are shown in light grey. For sake of clarity, they are only displayed for high divergence values (D_max>1µg m⁻² s⁻¹ or integrated D>0.1 kg s⁻¹), as they would otherwise dominate the figure for some regions like Europe, where 499 candidates are classified as “gap” and 1013 as “negative” (Table 3).

The map of the Middle East (a) contains most of the non-point-source examples shown in Fig. 3. Several more gaps occur all along the Persian Gulf coastline, and further negative candidates are found in mountains as well as along the coastlines.

Several large cities like (a) Cairo and Jeddah, (b) Paris and Madrid, (c) Delhi and Mumbai, or (d) Moscow and Saint Petersburg are categorized as area source. However, there are also some candidates categorized as area source which do not correspond to a megacity. In particular the candidate corresponding to the maximum divergence over India, which is caused by the coal-fired 5 GW Vindhyachal Super Thermal Power Station, was categorized as area source, as it interferes with the 4 GW Anpara power plant about 16 km northeast (Fig. 3d). Such sources could still be investigated in detail based on the divergence map. However, for interfering sources so close to each other, a quantification by a fully automated algorithm is challenging.

For Riyadh, the power plants PP9 and PP10 northeast and southeast of the city center are identified as point sources (compare Beirle et al., 2019). In contrast to Beirle et al. (2019), PP8 west of Riyadh is not identified as a point source as it could not be separated from Riyadh city, as a consequence of the strict pre-selection of candidates and the slightly larger fit interval, which were necessary in order to run the algorithm fully automated globally.

Several point sources are detected in the Middle East. There is a remarkable cluster of several point sources detected south of Baghdad. Note that there was even a point source detected within the Persian Gulf, which corresponds to the Zakum offshore oilfield.

The Indian subcontinent reveals the highest number of point sources of the investigated regions, contributing one-fourth of the global number. This reflects the quickly growing industrial activities, while measures for emission reduction still need to catch up. In addition, the divergence method obviously works very well for India, as the noise in divergence is quite low. This might be related to the dry season providing very good observation conditions without gaps, thereby suppressing sampling effects (compare Sect. 5.1.1).

In Europe, only very few point sources are detected, like the world's largest charcoal power plant Bełchatów in Poland; the German charcoal power plants Jänschwalde, Boxberg, and Neurath/Niederaußem (compare Beirle et al., 2019); or Europe's largest steel plant in Taranto, southern Italy. Remarkably, almost no point sources are detected for England and the Benelux countries, where the mean TVCD has a local maximum (Fig. 1a). Instead, there are several candidates classified as “gaps” and “negative”, which is related to the high noise observed in the divergence map. A similar situation is found for China, where TVCD is highest globally (Fig. 1a), but noise in D is large (Fig. S11 in the Supplement) such that the number of negative candidates is high (1013), but only few (34) point sources could be clearly identified. In Ukraine and western Russia, where mean TVCD levels are moderate, several point sources could be clearly identified. These striking regional differences in the performance of the automated point source detection will be discussed in detail in Sect. 5.1.

4.2 Catalog of point sources

We derive a global catalog of NO_x emissions from the 451 peaks categorized as point sources by sorting them according to the fitted emissions. A complete list of all detected point sources, including latitude/longitude and the estimated emissions, is provided in CSV format at https://doi.org/10.26050/WDCC/Quant_NOx_TROPOMI (Beirle et al., 2020) and also added as a data file to the Supplement.

Table 4 lists a selection of the identified point sources with some additional information on the respective NO_x source. It contains the top 10 emitters of the catalog. In addition, every 100th rank is included in order to illustrate conditions for lower divergence levels. Divergence maps for the same selection are shown in Fig. 5, where also external information on NO_x sources have been added. In the Supplement, respective divergence maps and tables of regional top emitters are shown for all considered regions.

Figure 5 Divergence maps for the selected point sources listed in Table 4. In addition to the location of the fitted Gaussian peak (grey), also some external information on GPPD power plants as well as industrial facilities and cities is added (green).

Most of the point sources listed in Table 4 can actually be associated with single or groups of power plants. Overall, a GPPD match was found for 242 point sources. The median distance between GPPD and point source locations was found to be 1.6 km, which is better than TROPOMI resolution. For the selection in Table 4, we did some additional inquiry and could identify the power plants Medupi (no. 3) and Presidente Vargas (no. 300), both missing in GPPD, as probable NO_x sources.

Table 4 Extract of the point source catalog for the top 10 and every 100th rank. Power plant capacity, fuel type, and facility names are added for matches to GPPD. The last two columns are not part of the catalog but have been added manually for the presented selection in order to provide information on the likely NO_x source where no (or insignificant) GPPD match has been found. Divergence maps for the same selection are displayed in Fig. 5. As discussed in detail in Sect. 5.2.6, the given emissions are biased low. Respective tables for regional top emitters are listed in the Supplement for all considered regions.

¹ Primary fuel of the GPPD match with highest capacity within 5 km.
² GPPD names have been shortened.
³ Missing in GPPD.

Other point sources are the Secunda CTL coal liquifying facility (no. 2), steel work facilities (no. 6, no. 10), and a cement plant (no. 100). Point source no. 7 is located in Ulsan, an industrial hotspot in South Korea. Here, however, we could not identify a single dominating NO_x point source. For the Hermosillo power plant, the peak fit is probably affected by Hermosillo city nearby (0.7 million inhabitants).

The four highest point source NO_x emissions are all found for South African coal power plants, which have already been presented in Beirle et al. (2019). Note that the emissions in Table 4 are lower than those given in Beirle et al. (2019) for various reasons, as discussed in detail in Sect. 5.2.6, mainly due to the missing AMF and lifetime corrections in the current study.

The thresholds for artifacts in divergence have been defined rather strictly. Consequently, the remaining locations listed in the catalog actually indicate stationary NO_x emissions. In the spot tests investigated exemplarily, we found no indication for false signals in the catalog.

4.3 Comparison to power plant database

Figure 6a provides a scatter plot of power plant capacity and point source emissions. Respective regional figures for all considered regions are provided in Fig. S15 in the Supplement.

Figure 6 Correlation between point source emissions and power plant capacity. (a) Scatter plot for all matches globally, with color coding primary fuel. (b) Scatter plot for coal-fired power plants in Australia, Europe, and South Africa.

Note that a perfect correlation between power plant capacity and NO_x emissions cannot be expected, as the emissions per capacity strongly depend on fuel type and technology and are particularly modified if emission control measures like selective catalytic reduction (SCR) are applied. High emissions for power plants with low capacity probably indicate other dominating NO_x sources nearby, like for Baotou (no. 6), where a 0.2 GW power plant was matching the point source location, but emissions are probably mainly caused by the metal smelting facilities. Low emissions from high-capacity power plants probably indicate the installation of SCR or a recent reduction in capacity.

High correlations can be observed for some regions like South Africa and Australia. This is probably indicating that the GPPD entries are reliable for these regions, the level of power plant technology is regionally similar, and the divergence method works well here. Figure 6b displays the scatter plots for coal-fired power plants for Australia, Europe, and South Africa. The slopes indicate clear regional differences in the emissions-per-capacity ratio.

5 Discussion

In this section we discuss the limitations as well as the potential of the presented point source catalog, give an outline on possible applications, and an outlook on improvements of the catalog in a future update.

5.1 Limitations

When interpreting the presented point source catalog, the following limitations have to be kept in mind.

5.1.1 Missing point sources

The catalog is incomplete, as point sources might be missing due to the following reasons:

1. Considered pixels. Only latitudes between 61^∘ N/S are considered, as for higher latitudes, the strict SZA threshold of 65^∘ would result in poor statistics. In addition, the criteria for defining the selection mask are quite strict in order to reduce the amount of data to be processed. There might thus be NO_x point sources not included in the selection mask. However, based on maps of the mean TVCD, we see no indication for strong point sources outside the considered area defined by M (compare Fig. 1a and b).

2. Gaps in input data. The mean divergence map reveals persistent gaps at some coastlines, in particular around the Persian Gulf. These gaps are caused by the cloud algorithm, as cloud retrievals are challenging for the transition of dark ocean to bright sand. The situation will improve for an updated cloud product which will be based on a ground albedo with higher spatial resolution.

3. Artifacts in divergence. In cases of inaccurate wind fields, as over mountainous terrain, as well as for systematic violation of the steady-state assumption, like systematic diurnal cycles of wind direction, the divergence map reveals artifacts, i.e., patterns of high negative values for D which cannot be explained by the loss term S. These cases are classified as “negative” and are thus missing in the catalog.

4. Noise in divergence. The noise level of the divergence map reveals large regional differences, causing respective differences in the performance of the peak fit algorithm. Noise levels are particularly high over regions with generally high TVCD levels, like eastern China or western Europe. This is caused by sampling effects. For high TVCD levels, also fluxes are generally high. As the daily flux maps have gaps due to cloud masking, the mean fluxes reveal “jumps”. This effect is probably intensified by the gridding by interpolation, where a gap (= missing value) in the input data results in gaps for a substantially larger area in the gridded data, as interpolation requires information from all surrounding pixels. Note that due to day-to-day changes of wind directions, a far higher amount of data would be needed in order to get smooth flux maps than to overcome the respective sampling issues in mean TVCDs. As the spatial derivative amplifies these jumps, the divergence is generally noisy over regions with high TVCD. Consequently, only few point sources are identified for the highly polluted regions in western Europe or eastern China, while many candidates are classified as “negative” due to the high noise levels.

5. Interfering sources. Point sources cause peaks in the divergence map which can be described by a Gaussian. Typical widths are ${\mathit{\sigma }}_{x,y}\approx \mathrm{5}$ km, in accordance to the spatial resolution of TROPOMI. Thus, point sources within a distance of less than about 20 km cannot be clearly separated by the automated algorithm. In cases of point sources close to each other, they are identified and quantified as one single point source (see next section), and their emissions are just added. If the distance is about 15–20 km, however, the joined peak from both sources is still processed as a single candidate but classified as “area source” due to the large width of the peak. For example, the Indian power plants Vindhyachal and Anpara located within 16 km (Fig. 3d) or PP8 at the western edge of Riyadh megacity (compare Fig. 2 in Beirle et al., 2019) are classified as “area source” and thus missing in the catalog.

5.1.2 Multiple sources

Due to the spatial resolution of TROPOMI, sources within about 10 km cannot be separated by the automatic peak-finding algorithm, and power plant stack emissions would automatically be merged with emissions from on-site activities such as heavy-duty diesel around the power plant. Several entries in the catalog contain multiple sources, like no. 1 (power plants Matla and Kriel within 5 km) or no. 3 (power plants Matimba and Medupi within 7 km). Similarly, the emissions from large industrial complexes, like the Ulsan industrial area, cannot be further specified. Thus the term “point source” means with respect to the TROPOMI spatial resolution. This is however still better resolved than many bottom-up inventories and typical horizontal resolutions of regional chemical transport models (CTMs). Thus, for emission inventories, also multiple and extended sources could appropriately be treated as point source emissions if used for models with a spatial resolution coarser than that of TROPOMI.

5.2 Uncertainties and accuracy

The main goal of v1.0 of the point source catalog is the identification rather than the quantification of NO_x point sources. Still we discuss and try to quantify the various sources for uncertainties of the derived point source emissions.

5.2.1 Gridding

Gridding is done by 2-D linear interpolation, which avoids abrupt jumps at the satellite pixel edges. Such jumps would cause large response in the divergence, i.e., the spatial derivative.

For narrow plumes, however, linear interpolation carries the risk of introducing a low bias. We have investigated this exemplarily for the power plant plume of PP9 northeast of Riyadh on 17 December 2017 (Beirle et al., 2019, Fig. 1a therein). The average plume TVCD from interpolation yields almost the same value as for conventional gridding (1 % lower). For the peak maximum, however, which is more relevant for the divergence than the mean, interpolation is 6 % lower. Thus, there is a low bias caused by linear interpolation, which is small compared to the other effects discussed below.

5.2.2 Lifetime

The quantification of NO_x emissions is based on the peak in D, ignoring the chemical loss of NO_x during downwind transport. We estimate the impact of neglecting the lifetime correction by comparing the catalog emissions to the respective emissions based on E (assuming a lifetime of 4 h). On average, the latter were higher than those based on D by about 25 %, which is quite small compared to other effects, in particular the AMF (see below).

In cases of much shorter lifetimes, like 1.5 h, as reported by Goldberg et al. (2019) for the Colstrip power plant in the USA (Fig. 1 therein), the emission estimates after lifetime correction would be higher by a factor of 1.73 on average.

5.2.3 Wind fields

As discussed in Beirle et al. (2019), different effects on wind field uncertainties have to be considered:

• Errors of the wind direction (both random and systematic) result in a underestimation of the flux and thus the estimated emissions, as any mismatch in wind direction leads to a low bias of the wind speed component projected to the actual wind direction. The underestimation is thus proportional to the cosine of the wind direction. For wind direction errors of 25^∘, it is about 10 %.

  Larger systematic errors in wind direction would cause visible artifacts in the divergence map and would thus be captured and removed by the check for negative divergence during candidate classification. In cases of larger random effects, the artifacts of individual days would at least partly cancel in the mean flux. But the wind speed component in the actual wind direction would be significantly underestimated, as well as the resulting divergence.

• The calculated fluxes, and thus the divergence map, are proportional to wind speed. The choice of the altitude of input wind fields thus affects the resulting emissions, as wind speeds are generally higher for higher altitudes. In Beirle et al. (2019), the effect for taking wind fields from 250 or 730 m rather than from 450 m on emissions was quantified as about 10 %.

  In this study, the focus was set to the clear identification of point sources. Thus, winds from a lower altitude (300 m) compared to Beirle et al. (2019) (450 m) were chosen in order to have optimal wind direction data close to the NO_x emission source. As discussed below, this choice removes the artifacts of divergence around the Matimba/Medupi power plants seen in Beirle et al. (2019).

5.2.4 Peak fit

The fitted emissions depend on the settings for fit radius and the forward model function. The fit window of 22 km around the source was chosen as compromise in order to cover the peak caused by a point source as far as possible while avoiding interference with neighboring sources. Variations of the fit settings within reasonable limits cause an uncertainty of about 20 % in the fitted emissions (Beirle et al., 2019).

5.2.5 AMF

Validation studies report on NO₂ TVCDs from TROPOMI being biased low for polluted sites (Verhoelst et al., 2021). This is probably due to the a priori vertical profiles, the cloud height product being biased low (Compernolle et al., 2021), and the surface albedo being biased high (Griffin et al., 2019), all resulting in high biased AMFs.

In Beirle et al. (2019), a correction of the low bias was applied by assuming the NO_x profile close to point sources to be in the lowest layer. Note that since the divergence is sensitive for the change of the NO_x flux due to the point source emissions, the required AMF correction has to be applied to the profile of the added rather than the total NO_x.

Still, we do not apply such a correction here, since there are indications that the cloud heights and surface albedo used for the NO₂ retrieval are systematically biased. Thus, also the provided AKs are biased, and a simple a posteriori correction of TVCDs is not possible.

For Germany, the AMF correction applied in Beirle et al. (2019) was a factor of 2 due to profile shape. In combination with a bias in the cloud height and ground albedo used for the TVCD retrieval, even higher factors have to be expected. The actual number will depend on surface albedo, aerosols, and cloud statistics (within the selection of low effective cloud fractions) and is high over dark surfaces and frequent cloud contamination but low over bright surfaces and few clouds (like for Riyadh). For partly clouded scenes, the AMF for the added NO₂ column could be easily too high by an order of magnitude if the real cloud is above the plume, but the retrieved cloud is below. Consequently, the added NO₂ would be biased low by an order of magnitude. This issue was improved in a recent update of the TROPOMI NO₂ L2 processor now using improved albedo and cloud products and thus more appropriate AKs. However, so far, no reprocessed NO₂ time series is available.

5.2.6 Low bias

As listed above, many effects contribute to the uncertainty of emission estimates from the divergence of mean NO₂ fluxes. Several of these effects are systematic in nature, resulting in an overall low bias of the derived emissions. In particular the effects of biased cloud height and inaccurate a priori profiles are expected to reveal large regional dependencies. Consequently, the low bias is hard to quantify for the global catalog. Thus we decided not to try to correct for the low bias of catalog emissions in this study but present the low biased estimates as they are with a clear disclaimer that the given emission estimates are biased low.

Accordingly, the emission estimates for South Africa reported in Beirle et al. (2019) were higher than the values listed in Table 4 by a factor of about 1.87 (for Matla/Kriel) up to 2.56 (for Medupi/Matimba). This discrepancy is a consequence of the missing AMF correction (factor 1.35 applied in Beirle et al., 2019), the missing lifetime correction (factor 1.1), the wind input data from lower altitudes (factor 1.1), and differences in grid definition, fit function (no rotation), and fit settings (factor 1.2), which together explain a factor of 1.9 as found for Matla/Kriel.

For Medupi/Matimba, the reason for the remaining discrepancy is the difference in input wind fields. The winds from higher altitudes (450 m) used in Beirle et al. (2019) cause a pattern of high negative D southwest of the power plants, indicating a mismatch of wind direction (see Fig. S2B in the Supplement of Beirle et al., 2019). Even in E, after applying lifetime correction, the negative values southwest remain in maps of E (Fig. 4A in Beirle et al., 2019). Consequently, the fit function finds a low background with a linear increase from southwest to northeast, and the emissions fitted on top of this low background are biased high.

In the current study, this artifact is not present (see point source no. 3 in Fig. 5), indicating that the wind direction from lower altitude (300 m) matches better to the actual NO_x transport of Medupi/Matimba emissions.

In the Supplement, catalog emissions are exemplarily compared to emissions reported by EPA for the top 10 emitters of the USA. A total of 7 of these 10 emitters are listed in the catalog, with correct naming found from the merging of GPPD. The other 3 emitters were also found as candidates but were classified as “negative”. EPA emissions were found to be higher than the emissions listed in the catalog by a factor of 3 (Navajo) up to 8 (Hunter). These power plants, however, are quite remote from large cities. Thus, in the absence of other significant NO_x sources, the modeled profiles used as a priori for the calculation of AMFs do not reflect the near-surface power plant plume due to the coarse model resolution. In addition, the cloud altitude used for the calculation of averaging kernels is biased low (Compernolle et al., 2021), causing high biased AMFs and low biased columns. The impact of this bias can easily reach 1 order of magnitude for cases where the retrieval assumes a cloud below the plume, while it is actually above. For these reasons, we have to expect that the low bias of the partial column added by a point source close to the ground is considerably larger than that of the complete tropospheric column, where validation studies report typical low biases of a factor of 2 for polluted sites. We will focus on this issue more quantitatively when preparing an update of the catalog, which we plan to process after reprocessing of the TROPOMI NO₂ product based on improved albedo and cloud products and thus more appropriate AKs.

5.3 Potential

Though the catalog is incomplete and the derived emissions are biased low, it still has the potential to improve our knowledge on NO_x emissions. Here we discuss the benefits of the divergence approach and the point source catalog and propose future applications.

5.3.1 Localization of NO_x and CO₂ sources

The Gaussian fit determines the location of the peak maximum. For all catalog locations which were inspected manually (i.e., all point sources listed in Table 4 or labeled in Fig. 4), Google Maps quickly reveals a plausible origin of the emissions close to the fitted peak location. For the examples shown in Fig. 5 which are related to power plants, the fitted point source location agrees to the actual power plant locations within about 2–3 km. Thus, the point source catalog accurately lists the location of NO_x point sources. As these sources are all related to combustion, this also provides valuable information on the location of CO₂ sources, which may also help to quantify CO₂ emissions from NO₂ measurements (Liu et al., 2020) and current and future satellite measurements of CO₂.

5.3.2 Spatial patterns

Within a regional focus, the catalog reflects the spatial distribution and relative importance of point sources. Regionally, high correlation between NO_x emissions and power plant capacity was observed (Fig. 6). Thus, it might be possible to re-distribute emissions from bottom-up inventories regionally according to the location of point sources in the catalog.

In addition, striking discrepancies between bottom-up inventories and the catalog, like strong point sources present in one but missing in the other, might be investigated in more detail and should result in improved emission inventories.

5.3.3 Up to dateness

Bottom-up inventories based on fuel consumption statistics have to collect and process input data from national reports. Thus, they have a time lag and are outdated when released for countries with high dynamics in industrial activities. Based on the divergence of NO₂ fluxes, single power plants can be identified and quantified, e.g., on an annual basis. This would allow us to detect short-term trends and power plants being switched on or off even in the vicinity of polluted cities such as Riyadh.

5.4 Outlook

This paper describes v1.0 of the NO_x point source catalog. We plan to update the catalog as soon as a reprocessed TROPOMI NO₂ product is available. We expect that some of the persistent gaps along coastlines, e.g., around Dubai, can be closed in the future, as soon as the cloud information is based on high-resolution maps of the surface albedo, and thus NO₂ TVCD becomes available there. In addition, improved surface albedo, cloud height, and a priori profiles are expected to improve the TVCD and (at least partly) remove the low bias.

In addition, we plan to use meteorological data from ERA-5 on higher spatial and temporal resolution. Currently, 6-hourly model output of ECMWF meteorological data was interpolated to a regular horizontal grid with a resolution of 1^∘. Case studies will be performed in order to find out which is the best compromise between spatiotemporal sampling and processing time.

6 Data availability

TROPOMI NO₂ data are available via Copernicus Sentinel-5P (processed by ESA), 2018, TROPOMI Level 2 Nitrogen Dioxide tropospheric column products, Version 01, European Space Agency, https://doi.org/10.5270/S5P-s4ljg54 (European Space Agency, 2019). The full NO_x point source catalog v1.0 is available at https://doi.org/10.26050/WDCC/Quant_NOx_TROPOMI (Beirle et al., 2020). The corresponding divergence maps are provided by the authors on request.

7 Conclusions

The high spatial resolution provided by TROPOMI allows for the detection and quantification of strong NO_x point sources like large power plants, metal smelters, cement plants, or industrial areas. We present v1.0 of a global catalog of NO_x point sources derived from the divergence of the mean NO₂ flux 2018–2019 by a fully automated iterative algorithm, yielding 451 point sources. A total of 242 of these point sources could be matched to combustion power plants from a global database, with a median distance of 1.6 km. The top four point source emissions are all located in South Africa and related to coal burning. About one-fourth of all point sources were found over the Indian subcontinent, where the method works quite well due to low noise levels.

The catalog is incomplete and misses point sources due to gaps in the divergence map (caused by gaps in the cloud product), artifacts in the divergence map (caused by non-steady state and inaccurate wind fields), noise in the divergence map (caused by sampling effects for regions with high background TVCD like western Europe or China), and interference of sources within about 10–20 km distance.

The listed emissions are biased low mainly due to a low bias of input TVCDs from TROPOMI. As this bias is expected to vary regionally and is hard to quantify, it is not corrected for v1.0 of the catalog. Exemplary comparisons to emissions reported from in situ monitoring reveals a low bias which can be as high as a factor of 8 for some power plants, which is probably caused by inappropriate a priori profiles and a low bias in the cloud height.

Still, the catalog has high potential for checking and improving emission inventories, as it localizes NO_x (and also CO₂) sources, provides spatial patterns of the distribution of sources, and yields up-to-date emission data.