Catalog of NOx emissions from point sources as derived from the divergence of the NO2 flux for TROPOMI

We present version 1.0 of a global catalog of NOx emissions from point sources, derived from TROPOMI measurements of tropospheric NO2 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 NOx point source by a fully automated algorithm, while 5 ambiguous cases as well as area sources such as Megacities are skipped. 242 of these point sources could be automatically matched to power plants. Other NOx point sources listed in the catalog are metal smelters, cement plants, or industrial areas. The four largest localized NOx emitters are all coal combustion plants in South Africa. About 1/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 NO2 product at some coastlines, inaccurate or 10 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 NO2 in TROPOMI NO2 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 NOx, and thus also CO2. 15 The catalog of NOx emissions from point sources is freely available at https://doi.org/10.26050/WDCC/Quant_NOx_ TROPOMI (Beirle et al., 2020). Copyright statement. This work is distributed under the Creative Commons Attribution 4.0 License.

spatial resolution are essential (Bouarar et al., 2019). Such data is 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 2 and other 25 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 2 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 2 profile shape as well as on viewing geometry, surface albedo, aerosols, and particularly on clouds. 30 NO 2 TVCDs from satellite measurements provide independent information on the spatial distribution and strength of tropospheric NO 2 levels on global scale since the mid nineties, 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 single 35 payload of ESA's Sentinel-5 Precursor satellite mission. TROPOMI provides NO 2 TVCDs on unprecedented high spatial resolution (7.2×3.6 km 2 until 5 August 2019, 5.6×3.6 km 2 thereafter) and with a high signal to noise ratio (Geffen et al., 2020). Single TROPOMI overpasses clearly reveal NO 2 plumes downwind from strong NO x sources like large power plants . In temporal mean NO 2 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). Upscaling NO 2 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: 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 .
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 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 60 TVCDs.
From the derived divergence maps, a catalog of NO x point sources is extracted by a fully automated algorithm.
The manuscript is organized as follows: In section 2, the input data sets used in this study are specified. The detailed data processing is explained in section 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.

65
2 Input data 2.1 Tropospheric NO 2 column densities The NO x point source catalog is based on NO 2 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 copernicus.eu. Details of the TROPOMI tropospheric NO 2 product are given in Geffen et al. (2019) and Table 1. Processing settings in this study as compared to Beirle et al. (2019).

Ozone climatology
Ozone mixing ratios, used for the scaling of NO 2 to NO x , were taken from the Earth System Chemistry integrated Modelling 85 (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 OMI-Aura overpass time (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 NO x /NO 2 ratio is that close to ground. We thus took O 3 concentrations from the lowest model layer.

Power plant database
The Worlds Resources Institute provides an open access Global Power Plant Database (GPPD) (Byers et al., 2019). We use this database 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 95 and oil as primary fuel, respectively.

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). For this study, we select TROPOMI tropospheric NO 2 column densities V NO2 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 SZA above 65 • for the calculation of fluxes. This strict criterion removes obser-105 vations for sun being low, and implicitly results in a gradual removal of wintertime measurements for mid-latitudes, 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 to simply parameterize the NO 2 photolysis as function of the SZA (see section 3.4).
Since large parts of the globe are free from stationary NO x point sources, in particular oceans, deserts and forests, the processing focusses 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 2 TVCDs. The constrution of the selection mask is explained in detail in the Supplement.

115
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 NO2 are 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 2D interpolation of TROPOMI pixel centers using the griddata function from the 120 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.

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 the focus of this study are emissions from point sources, we consider wind fields representative for the transport of 130 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 minutes at 60°S and 60°N at solstice ± 7 minutes for summer/winter. Thus the orbit time stamp reflects the wind conditions for recent plume histories.
135 Table 2. Definition of regions used for the regional statistics shown in table 3 and for regional figures shown in the Supplement. 3. in lat/lon to the 0.025°grid.

Up-scaling of NO 2 to NO x
In this study we scale the measured NO 2 TVCD to a NO x TVCD for each TROPOMI pixel. The conversion factor L is calculated according to the photostationary steady-state

140
The impact of 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 2 J is parameterized as function of the SZA θ by 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).

145
The reaction rate constant k for the reaction of NO with O 3 is parameterized as function of temperature (in Kelvin) by k = 2.07 × 10 −12 × exp(−1400/T ), the NO is converted to NO 2 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 (section 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.

Gridded fluxes and divergence
From gridded NO x columns and gridded wind fields, the gridded NO x flux in both x and y direction is derived for each 165 TROPOMI orbit. Mean fluxes are calculated for the period 2018-2019, where calm wind conditions (w < 2m/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 implicitely by the strict SZA threshold of 65 • .
From the mean zonal and meridional flux maps, the divergence map D = ∇ · F is calculated, which is the basis for the 170 identification and quantification of point sources below.

Lifetime and background corrections
In Beirle et al. (2019), emission maps were derived by adding the sink term S = V /τ to the divergence map. For this, a constant lifetime of τ = 4 h was applied, as derived from OMI data for Riyadh . In addition, V was corrected by subtracting the regional background, as the lifetime estimate was derived for freshly released, surface near pollution, while 175 upper tropospheric background NO 2 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 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:

180
-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, volantile organic compounds (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 185 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).

190
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 /τ 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.

AMF correction
Tropospheric column densities of NO 2 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 sum of height-dependend "box AMFs", representing the vertical measurement sensitivity, weighted by the relative profile (Wagner et al., 2007). In the operational TROPOMI NO 2 product, averaging Kernels (AKs) are provided that are proportional to box-AMFs and allow to correct the AMF for a 200 different vertical profile (Eskes et al., 2003).
Validation studies report on a general low bias of NO 2 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  and the albedo maps used are biased high (Griffin et al.,205 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 surface-near 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 2 product.
Consequently, the low bias of TROPOMI TVCDs is directly transferred into the NO x emissions listed in the catalog, as 210 discussed in detail in sect. 5.2. The low bias is expected to be improved with an updated NO 2 product, which will then be used for deriving an updated version of the point source catalog.

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 215 skipped. Thus, the resulting catalog of point sources is incomplete; a detailed discussion on various reasons for missing point sources is given in section 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 as point source candidate. 220 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 −2 s −1 , which is < 4% of the initial 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 to decrease this threshold and investigate additional small NO x point sources.

Pre-classification of candidates 230
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 zooms of the divergence map 30 km around the candidate. As soon as the candidate is classified in one of the following categories, the pre-classification stops. In Fig. 3, examples for each category are shown.

Category gap:
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 3(a) displays 240 an example for a candidate of category "gap" around Dubai, where the global maximum of divergence was found, but missing data does not allow for further quantifications.

Category negative:
If the minimum divergence within the 30 km zoom is negative with an absolute value larger than 50% of the maximum value, the candidate is "negative".

245
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 re-moval is taking place over large distances (at wind speed of 5 m/s, a lifetime of 4 hours 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 250 lifetimes), but indicate an inappropriate simplification of the complex 3D wind fields by a 2D 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 3(b) 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 2 around Tehran that can even be spotted in the global map 255 ( Fig. 1(a)). 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.

Category area source:
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 half

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: with the fit parameters A being the peak integral, σ x , σ y the width of the Gaussian, ∆x, ∆y 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 275 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 fit := 2.355 × σ x,y exceeds 17 km for x or y.

280
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 detection limit for optimal conditions) or a relative fit error above 30% are categorized as "uncertain".

Candidate removal 285
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.

290
Depending on the classification, the following procedure is applied: -For point sources, an ellipse with 2×σ x,y as 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.

295
-For area sources, also 30 km around the maximum are removed, if classified based on W first guess . If the area source was classified based on W fit , an ellipse with 2×σ x,y is removed as for point sources.

Identification of point sources associated with power plants
We   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. 3 d). Such sources could still be investigated in detail based on the divergence 320 map. However, for interfering sources such 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 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.

325
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.  (Fig. 1 a). 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. 1 a), 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 West Russia, where mean TVCD 340 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.

Catalog of point sources
We  Table 4 lists a selection of the identified point sources with some additional information on the respective NO x source. It contains the top ten 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 350 have been added. In the Supplement, respective zooms of the divergence maps and tables of regional top emitters, are shown for all considered regions.
Most of the point sources listed in table 4 can actually be associated to 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 355 power plants Medupi (#3) and Presidente Vargas (#300), both missing in GPPD, as probable NO x source.
Other point sources are the Secunda CTL coal liquifying facility (#2), steel work facilities (#6, #10), and a cement plant (#100). Point source #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 M inhabitants).

360
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 section 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 365 false signals in the catalog.

375
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. Fig. 6 (b) 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.

380
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. Table 4. Extract of the point source catalog for the top ten 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 NOx source where no (or insignificant) GPPD match has been found. Zooms of the divergence map for the same selection are displayed in Fig. 5. As discussed in detail in section 5.2.6, the given emissions are biased low. Respective tables for

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

Missing point sources 385
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.   1(a) and (b)).

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.

395
The situation will improve for an updated cloud product which will be based on a ground albedo with higher spatial resolution.

Artifacts in divergence:
In case 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 400 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.

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 405 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 410 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.

Interfering sources:
415 Point sources cause peaks in the divergence map which can be described by a Gaussian. Typical widths are σ x, y ≈ 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 case 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 420 source due to the large width of the peak. For example, the Indian power plants Vindhyachal and Anpara located within 16 km (Fig. 3 d) 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.

Multiple sources
Due to the spatial resolution of TROPOMI, sources within about 10 km cannot be separated by the automatic peak finding 425 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 #1 (power plants Matla and Kriel within 5 km) or #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 430 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.

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 435 we discuss and try to quantifiy the various sources for uncertainties of the derived point source emissions:

Gridding
Gridding is done by 2D 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 exem-440 plarily 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 than 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 however small compared to the other effects discussed below.

445
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 case of much shorter lifetimes, like 1.5 h, as reported by Goldberg et al. (2019) for the Colstrip power plant in the USA 450 ( Fig. 1 therein), the emission estimates after lifetime correction would be higher by a factor of 1.73 on average.

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

475
Validation studies report on NO 2 TVCDs from TROPOMI being biased low for polluted sites . This is probably due to the a-priori vertical profiles, the cloud height product being biased low , and the surface albedo being biased high , 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 being completely in the lowest layer. Note that since the divergence is sensitive for the change of the NO x flux due to the point source 480 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 2 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 485 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 2 column could be easily too high by an order of magnitudes, if the real cloud is above the plume, but the retrieved cloud is below. Consequently, the added NO 2 would be biased low by an order of 490 magnitude. This issue was improved in a recent update of the TROPOMI NO 2 L2 processor now using improved albedo and cloud products, and thus more appropriate AKs. However, so far, no reprocessed NO 2 timeseries is available.

Low bias
As listed above, many effects contribute to the uncertainty of emission estimates from the divergence of mean NO 2 fluxes.
Several of these effects are systematic in nature, resulting in an overall low bias of the derived emissions. In particular the 495 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 500  (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

510
In the current study, this artifact is not present (see point source #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. 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 515 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 absence of other significant NO x sources, the modelled 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 , causing high biased AMFs and low biased columns. The impact of this bias can easily reach one order of magnitude for cases where the retrieval assumes a cloud below 520 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 2 product based on improved albedo and cloud products and thus more appropriate AKs.

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: CO 2 sources, which may also help to quantify CO 2 emissions from NO 2 measurements (Liu Fei et al., 2020) and current and future satellite measurements of CO 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 540 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.
Bottom-up inventories based on fuel consumption statistics have to collect and process input data from national reports. Thus, 545 they have a time lag and are outdated when released for countries with high dynamics in industrial activities. Based on the divergence of NO 2 fluxes, single power plants can be identified and quantified e.g. on annual basis. This would allow to detect short-term trends and power plants being switched on or off even in the vicinity of polluted cities such as Riyadh.

Outlook
This manuscript describes v1.0 of the NO x point source catalog. We plan to update the catalog as soon as a reprocessed 550 TROPOMI NO 2 product is available. We expect that some of the persistent gaps along coastlines, e.g. around Dubai, can be closed in future, as soon as the cloud information is based on high resolution maps of the surface albedo, and thus NO 2 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 555 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 spatio-temporal sampling and processing time.

Conclusions
The high spatial resolution provided by TROPOMI allows for the detection and quantification of strong NO x "point sources" The catalog is incomplete and misses point sources due to gaps in the divergence map (caused by gaps in the cloud product), 565 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 570 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 2 ) sources, provides spatial patterns of the distribution of sources, and yields up-to-date emission data. The full NO x point source catalog v1.0 is available at https://doi.org/10.26050/WDCC/Quant_NOx_TROPOMI . The corresponding divergence maps are provided by the authors on request.
Author contributions. SB performed the analysis and wrote the manuscript with input from all co-authors. CB processed the TROPOMI L2 data. SD processed ECMWF meteorological data. HE performed the retrieval of operational TROPOMI L2 NO2 TVCDs. VK compiled the Ozone climatology from ESCiMo model data. AdL initiated the matching between catalog point sources and GPPD. TW contributed to the 580 data analysis and interpretation and supervised the study.

Competing interests. None
Medium-Range Weather Forecasts (ECMWF). Andrea Pozzer is acknowledged for providing the ESCiMo model data used for the Ozone climatology.