Quantifying urban, industrial, and background changes in NO 2 during the COVID-19 lockdown period based on TROPOMI satellite observations

Abstract. The COVID-19 lockdown had a large impact on anthropogenic emissions of air pollutants and particularly on nitrogen dioxide (NO2). While the overall NO2 decline over some large cities is well-established, its quantification remains a challenge because of a variety of sources of NO2. In this study, a new method of isolation of three components: background NO2, NO2 from urban sources, and from industrial point sources is applied to estimate the COVID-19 lockdown impact on each of them. The approach is based on fitting satellite data by a statistical model with empirical plume dispersion functions driven by the observed winds. Population density and surface elevation data as well as coordinates of industrial sources were used in the analysis. The NO2 vertical column density (VCD) values measured by Tropospheric Monitoring Instrument (TROPOMI) on board Sentinel‐5 Precursor over 263 urban areas for the period from March 16 to June 15, 2020, were compared with the average VCD values for the same period in 2018 and 2019. While background NO2 component remained almost unchanged, the urban NO2 component declined by 18–28 % over most regions. India, South America, and a part of Europe (particularly, Italy, France, and Spain) demonstrated a 40–50 % urban emissions decline. In contrast, decline over urban area in China, where the lockdown was over during the analyzed period, was only 3 % except for Wuhan, where more than 60 % decline was observed. Emissions from large industrial sources in the analyzed urban areas varies largely from region to region from +5 % for China to −40 % for India. Changes in urban emissions are correlated with changes in Google mobility data (the correlation coefficient is 0.66) confirming that changes in traffic was one of the key elements in decline of urban NO2 emissions. No correlation was found between changes in background NO2 and Google mobility data.


NO2 VCD spatial distribution used for comparisons with satellite-based estimates. The analysis was done in three stages. The algorithms were originally developed and tested for the U.S. and Canada, then it was applied for Europe. For these regions, detailed information about emission sources was available. Finally, it was applied worldwide, where information about the emission sources was often limited.  where α0, αe, αp, αi, β0, β1, and β2 are the unknown regression parameters representing population density-related proxies and emissions from individual point sources and a background with contribution from the elevation; Ωp is the source plume function for the population density-related distributed source (or area source); Ωi are the source plume functions for industrial point sources; H is the elevation above sea level and the scaling factor H0=1 km was introduced to make the exponential argument dimensionless; and ε is the residual noise. The fitting was done for 3° by 4° areas around large cities by minimization of the 5 squares of the residuals (ε). The fitting and parameter estimation was done using all individual level 2 pixels for the period from March 16 to June 15 three times: for 2018, 2019, and 2020. So, three sets of coefficients were obtained and then used to estimate the background levels and emissions. Then, the average results for 2018 and 2019 were compared with these for 2020.
As in Fioletov et al., 2017, the plume from a point source i is described by a plume function Ω (θ, φ, ω, s, θi, φi) where θ and φ are the satellite pixel coordinates; ω and s are the wind direction and speed for that pixel; and θi and φi are the source 10 coordinates. An unknown parameter (αi) represents the total NO2 mass emitted from the source i. The emission rate for source i can be expressed as Ei =αi/τ, where τ is a prescribed NO2 lifetime (or, more accurately, decay time, but we use the term "lifetime" because it is more common). Note that τ is different from the chemical lifetime (de Foy et al., 2015). Once the emission rate is established, it can be used to reconstruct how distribution of NO2 emitted by that source would be seen by a satellite. 15 The plume functions Ω are EMG functions that are commonly used to approximate plumes of VCDs of trace gases such as NO2, SO2, and ammonia (Beirle et al., 2011(Beirle et al., , 2014Dammers et al., 2019;Fioletov et al., 2017Fioletov et al., , 2015de Foy et al., 2015;Liu et al., 2016;McLinden et al., 2020). Similar in concept to a Gaussian plume function, they also take into account the finite physical size of the source and the spatial resolution of the satellite instrument being utilized. The shape of the EMG function depends on the prescribed plume width (w) and the gas lifetime () plus on an unknown emission strength of the 20 plume αp. In general, the plume width parameter depends on the size of the source and the size of satellite pixel. The value of 8 km for plume width was used in this study for TROPOMI along with a constant lifetime value of 3.3 hours. The switch from 7 to 5.6 km along-track resolution in 2019 might have some impact on the optimal plume width, but as the sensitivity analysis show that small changes in w have only a minor impact on the results (Supplement A). The lifetime reflects the rate at which NO2 is removed from the plume due to chemical conversion or physical removal such as deposition; it depends on several 25 factors such as season and NO2 concentration. It is about 2-6 hours in summer and longer in winter (de Foy et al., 2014;Liu et al., 2016). Moreover, for some sources, the lifetime may be changing over time (Laughner and Cohen, 2019) as NO2 concentration declines, although other studies suggest that such changes are minor (Stavrakou et al., 2020). Recent TROPOMIbased estimates show that a typical lifetime in urban areas is between 2 and 5 hours in spring and autumn with shorter lifetimes at low latitudes (Lange et al., 2021). While the lifetime has a large impact on the emission estimates, relative changes are less 30 sensitive to it. The impact of the plume widths and lifetime parameters on the estimates and standard deviation of the residuals is discussed in the Supplement A.
Unlike many previous studies (Beirle et al., 2011;Fioletov et al., 2016;Lange et al., 2021) where the background offset was presumed to be constant and estimated from, for example, upwind NO2 data, we included a special term that is https://doi.org /10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. responsible for it. In equation (1), the α0 + (β0 + β1(θ -θ0) + β2(φ -φ0))·exp(-H/H0) term represents the "background" NO2 that is assumed to be declining exponentially with elevation, i.e., within the analyzed 3° by 4° area, the higher is the elevation the lower is tropospheric column NO2. It was also assumed that this contribution from elevation depends on geographical coordinates only and not on the winds. Even in absence of any sources, there could be some gradient in tropospheric NO2 over the analyzed area, as for example, over some regions in northern Canada or along the east coast of the U.S. (Fig. 1). To account 5 for such gradients, the linear term β1(θ-θ0)+β2(φ-φ0), where θ0 and φ0 are the coordinates of the centre of the analyzed area, was added. In other words, it was assumed that there is a linear gradient of background NO2 within the analyzed area and NO2 VCD declines exponentially with height over elevated regions. Finally, α0 was added to the model to account for remaining free-tropospheric NO2 at high elevations where exp(-H/H0) is very close to 0. Its presence gives a better agreement with the satellite data for areas with a high range of elevations. Since this term is a part of the statistical model, all parameters α0, β0, 10 β1, and β2 are estimated from the fitting. Once they are estimated, the term can be calculated for each location within the analysed 3⁰×4⁰ area that gives a "background" value for that location that depends on the coordinates and elevation only. For simplicity, we will refer to the term discussed in this paragraph as to the "background" component.
The ap Ωp term represents the emissions contribution from factors, related to human activity, where the population density is used as a proxy. The population density data were converted to a 0.2° by 0.2° grid and then each of 336 (16 x 21) 15 grid cells were considered as point sources located at the grid centre with emissions proportional to the population. However, instead of estimating emission for each such "source" separately, the emissions were parameterized by a single unknown parameter αp. This means that emissions per capita are assumed to be the same everywhere within the analysed 3° by 4° area.
The composite plume function Ωp is, therefore, a sum of plume functions of all individual 0.2° by 0.2° grid cell centres multiplied by the grid cell population. Thus, Ωp depends on geographical coordinates, population density and local winds. 20 Eq.1 is a linear regression statistical model with unknown coefficient sets α and β. The proxy functions used in the model preferably should be uncorrelated, because otherwise the coefficients have correlated errors making their interpretation difficult. There are three main components in the model: related to background and elevation (four fitted coefficients), to the population density (one coefficient) and to industrial sources (variable number of coefficients from zero to few dozens). We will refer to them as to background, urban and industrial components. For a typical urban area, these components should be 25 independent: high population density zones typically occupy a small part of the area and industrial sources are typically located away from such highly populated zones. Note that the lifetime is relatively short (3.3 hours) and the median wind speed in, for example, the eastern U.S. is about 10 km per hour, so sources located 30-40 km apart typically have uncorrelated plume functions.
High correlation between the population and landscape-related proxies is possible if a city is in a valley surrounded 30 by high mountains. The correlation could be reduced by increasing the size of the analysed area, but if the area is too large, the assumption that the background level has a linear gradient in the area may not be valid. Therefore, we limited the area to 3° by 4°. The correlation coefficients between the site elevation and population density for 3° by 4° areas are typically small.
For example, in the U.S., correlations are positive over Florida (about 0.2) with the population density higher in the area above https://doi.org /10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. sea level, and negative in the Portland-Seattle-Vancouver area (about -0.35), where it is higher near the ocean and lower in the mountains. As the plume functions of individual industrial sources are very local (~50 km footprint), they do not correlate with the elevation. With such low correlation coefficients, elevation does not affect estimates of other parameters of the regression model.
When industrial point sources are located in close proximity, their plume functions in the statistical model (Eq.1) are 5 highly correlated. In practice, it often appears as if, for example, estimated emissions from one source are unrealistically high, while emissions from the other near-by source are low or even negative. In such cases, emissions from individual industrial sources often cannot be estimated. However, the sources can be grouped into independent clusters and total emission from such clusters can be estimated. Such grouping could be done manually on a case-by-case basis, but it would be subjective and very time consuming. Instead, we applied an algorithm based on factor analysis. We would like to emphasize, that the factor 10 analysis, described in the next three paragraphs, was used to improve emission estimation for individual sources or cluster of sources. It is not required if only total emissions from all point sources in the area are estimated in order to separate them from urban emissions or if all industrial sources are isolated remote sources.
To group industrial sources into clusters, an orthogonalization process was applied to the plume function of individual industrial sources. First, the correlation matrix for the plume functions of individual point sources (Ωi) was calculated and 15 eigenvalues and eigenvectors (or "factors") of the correlation matrix were determined. The correlation matrix was calculated just once using March 16 -June 15 data from all three years. An isolated remote source would appear as an eigenvector with an eigenvalue of 1. Two (or more) closely located, but isolated from others, sources, would have one corresponding eigenvector and an eigenvalue of 2 (or more). Eigenvalues lower than 1 mean that the corresponding sources are already partially included in other eigenvectors. To reduce the number of "factors", only "factors" with eigenvalues > 0.6 were kept. 20 The approach based on eigenvalues of the correlation matrix creates proxies that are not correlated and reduces the number of the fitting coefficients. While they correctly describe the total contribution of all industrial sources in the area in the total NO2 variability (or total emissions), individual eigenvalues, i.e., linear combinations of the original plume functions, may not have clear interpretation. For example, they may include the original plume functions with negative coefficients. In order to avoid that and obtain proxies that have a meaningful interpretation, the eigenvalues were linearly transformed, so they 25 became as close to the original plume functions as possible, while the correlation coefficients between them remained low.
This was done using the varimax factor analysis method that is implemented in modern statistical software packages such as R and SAS (Belhekar, 2013). It orthogonally rotates the established "factors" to maximize the sum of squared correlations between the original variables and factors. Then, the algorithm uses a linear combination of the original variables that have the highest correlations with the rotated "factors". I.e., the condition of orthogonality is removed in order to find the 30 simplest linear combination of the original variables. In practice, the algorithm produces a set of "clusters", i.e., linear combinations of the original plume functions, that have low correlation coefficients (typically less than 0.2) between them, and each cluster has high correlation coefficient (typically more than 0.95) with one orthogonal "factor". To simplify this further, if a linear combination has a weight for an original variable under 0.2, its weight was set to 0. As a result, all non-https://doi.org/10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. isolated point sources were grouped into small clusters and emissions estimates were done for such clusters instead of individual sources, while each isolated remote source forms a single-source cluster that corresponds to only that source. It is possible that a single source contributes to more than one cluster that makes interpretation of emissions for such clusters more difficult, but such cases are rare.
As in any regression analysis-based study, correlation between the proxies is one of the main obstacles in the result 5 interpretation. The "orthogonalization" of plume functions from industrial sources largely reduces cross-correlations between the proxies, but high correlations between industrial and population density-related plume functions are still possible if industrial sources are located in highly populated areas. In such cases, it may be difficult to separate its signal from the contribution of the population density-related proxy. For example, in one case (Edmonton, Canada) this correlation coefficient was as high as 0.94 and it was not possible to separate urban and industrial emissions there. For seven cities in the U.S. and 10 Canada it was between 0.81 and 0.85 (for the remaining sources, it is less than 0.77). Note that for large cities and small industrial sources, high correlation means that the industrial emissions cannot be reliably estimated, although the impact on estimation of the population density-relate signal is small. We monitored the correlation coefficients between industrial and population density-related plume functions and, in some cases, excluded certain sources or even certain urban areas from the analysis. 15 The algorithm is illustrated in Fig. 2, where individual terms of Eq. 1 are shown for an area centred on Montreal. The area includes two large cities, Montreal (4.2 million) and Ottawa (1.4 million, including the sister city of Gatineau). The terrain elevations in the analysed area are in the range from just a few metres above sea level along the Saint Lawrence River to more than 500 meters 100 km north of Montreal. For this plot (as well as for Fig. 1 and other figures), we used a non-linear scale that is more sensitive to small quantities in order to make small deviations more pronounced. The top row of Fig. 2 shows the 20 mean TROPOMI NO2 data (Fig. 2a), the fitting results (Fig. 2b), and the difference between them or the residuals (Fig. 2c).
Individual fitting components such as the background term (α0 + (β0 + β1(θ -θ0) + β2(φ -φ0))·exp(-H/H0)), population density- The suggested algorithm essentially finds the emission levels that give the best agreement with the TROPOMI data and then uses these estimates to "reconstruct" the spatial NO2 distribution as well as contribution from each source. As explained by (Fioletov et al., 2017), the technique of VCD reconstruction from fitted coefficients αi using Eq.1 used to isolate different components, can be applied to the reported emissions Ei. by using αi= Ei×τ. This produces a map of VCD that would be seen by satellites if these reported emissions are the only sources of NO2. The same approach was employed here. In the 20 U.S., directly measured industrial NOx emissions from point sources are available. It was assumed that the NO2 to NOx ratio is constant 0.7. The map of NO2 from the reported emissions is shown in Fig. 4 column (g). We would like to emphasise that such reconstruction is based on emissions data only, without any satellite NO2 observations (although τ and w in the plume functions were the same as in the satellite-based estimates). Finally, the maps of the difference between industrial sourcesrelated component and NO2 from reported emission-based reconstruction is also shown (column h). 25

Case studies
In the case of Boston, there is a single urban source with no large industrial sources nearby and with relatively small impact from the terrain. The Atlanta area represents the case where the urban component is well-separated from industrial sources. In the Pittsburgh area, industrial and urban sources have comparable contributions, and the industrial emissions estimates can be validated by measured emissions. Multiple industrial sources in the Huston area are missing from the used EPA NEI emission 30 https://doi.org/10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. database and this example illustrates why it was necessary to add locations of additional sources to explain the observed NO2 distribution in the area.
Boston is a major urban area with a population of more than 8 million (for the Combined statistical area of Greater Boston). On the TROPOMI NO2 map (Fig,4, col. a), it appears as a large "hotspot" that can be successfully reproduced by the statistical model (Eq.1) using the population density as a proxy. There is a 25% decline in the urban emissions in 2020 5 compared to the 2018-2019 average. It also shows one of the largest in the U.S. declines in the background component (about 15%).
The Atlanta area hosts the Hartsfield-Jackson Atlanta International Airport, the world's busiest with more than 100 million passengers per year in 2018-2019 (https://aci.aero/data-centre/annual-traffic-data/passengers/2017-passengersummary-annual-traffic-data/). The Atlanta airport NO2 signal can be easily isolated since the airport is located far away from 10 industrial sources (the correlation coefficients between the plume functions are less than 0.2) and on a distance from Atlanta city's most populated area (the correlation coefficient is 0.54). VCDs estimated for the industrial source clusters (column f in  Fig.4). When examining their difference (column h), the airport appears as a hotspot as its emissions are not included in the reported emissions inventory. All industrial sources (power plants) are also far away from the populated area and signals from these sources can be clearly separated from the 15 population density-related signal. The area is flat, and the background level does not change very much within the area except for mountains to the north.
The Pittsburgh area has one of the highest emissions from industrial sources among all analysed areas in the U.S.
Several coal-burning power plants are located east, west, and south of the city. The emissions from them are comparable or even larger than from the city itself. As mentioned, the NO2 distribution around major industrial sources can be "reconstructed" 20 from the emissions, reported by these sources (column g in Fig.4). The main features of the reconstructed NO2 distribution from industrial sources based on satellite estimates, the reported emissions are rather similar, and the differences (column f minus column g) are small, although NO2 from reported emissions are slightly larger for the cluster of power plans east of the city. The background values are fairly constant within the area except for the Appalachian Mountains area.
In the case of Houston, the EPA NEI emissions inventory contains emissions from the power plants in the area, but 25 not from large oil refineries that are responsible for hotspots seen on the TROPOMI mean NO2 plot. In order to reduce the residuals, we added the coordinates of 10 known refineries as industrial sources in Eq. 1. As a result, the residuals became small, although some elevated NO2 values can still be seen in the industrial areas suggesting the list of refineries included in this analysis may not be complete. This example illustrates that it is necessary to add some sources to the list of sources based on the used EPA inventory, derived from measured emissions. A total of 56 sources across the US and Canada were added 30 and the list of sources in the Supplement (file Additional_Sources_Canada_US.xls).
Landscape elevation plays an important role in the NO2 distribution that can be illustrated by the following examples ( Fig. 5). If the surface is nearly flat in the analysed area (as, for example, in the case of Minneapolis), the background component is dominated by a linear gradient. However, the elevation affects the NO2 distribution near mountain areas as, for https://doi.org/10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. example, in case of Seattle where mountains as high as 2000 m are located east of the city. It is interesting to note that the background components are practically identical for both periods that gives a high confidence in the obtained results. The influence of the landscape on the NO2 distribution also explains why the distribution near Seattle does not look like a "hotspot" NO2 distribution near a typical large urban area.

Relative contribution of different components 5
NO2 VCD represents total number of molecules, and equivalently mass, per area unit. When individual components affecting the NO2 distribution are known, it is possible to estimate their relative contribution to the total mass based on Eq.1. The diagram in Fig. 6 shows such contribution of individual components for the Montreal area. Most NO2 mass is associated with the term related to the background component. For the Montreal area, the contribution of industrial sources is four times less than the contribution of the urban component and the two of these components are responsible for less than a quarter of the 10 total NO2 mass in the area.
The mean NO2 distribution near major emission sources has sharp gradients and estimates of the NO2 lifetime from satellite data suggest that this time is relatively short, on the order of a few hours (Beirle et al., 2011;de Foy et al., 2015). It is even shorter than for satellite SO2 VCDs (Fioletov et al., 2015). However, large background component may suggest that the lifetime should be relatively long since NO2 distribution follows the terrain over large areas. This difference in the lifetime 15 could be reconciled if we assume that a fraction of NO2 emitted from cities and industrial sources gets into free troposphere and have a longer lifetime there than near the ground. Also, levels of the OH radical, the main chemical NOx sink, within a plume can be much larger than under "clean" conditions and NO2 lifetime could be longer under such condition than in the plume (Juncosa Calahorrano et al., 2021). Other sources, e.g., lightning or soil emissions may contribute to background component NO2 directly. The background term can also include components of stratospheric NO2 that was imperfectly 20 removed as part of the retrieval algorithm.
Relative contribution of the three components for the 2018-2019 period are shown in Fig, 7. Plumes from urban and industrial sources, or to be precise, αp Ωp and Σ αiΩi terms in Eq.1, are responsible for less than a third of total NO2 mass in all of the analysed 3° by 4° urban areas in the U.S and Canada. Most of NO2 mass belongs to the background component that is not directly linked to particular urban or industrial sources. Fig. 7 also shows that NO2 mass emitted from cities are larger than 25 emissions from the industrial sources for most of analysed areas in the U.S. and Canada.

The COVID-19 lockdown impact: the U.S. and Canada
The ability of the method to isolate individual components of the total NO2 mass makes it possible to estimate the impact of the COVID-19-related lockdown on these components. The results for the background, urban, and industrial components are shown in Fig. 8 To illustrate the changes in the background component, Fig. 8 (top) shows the mean VCD values of that component shown in Fig 4, column d (or, in other words, the mean value of α0 + (β0 + β1(θ -θ0) + β2(φ -φ0))·exp(-H/H0)) for the analyzed areas for the two time intervals (left) as well as the percentage change in 2020 vs. 2018-2019 values (right). The mean value of decline for the background components among all urban areas is -1.0% ± 3.3% (all error bars in the study correspond to 2σlevel). Thus, the average background value in 2020 was almost the same as in 2018-2019. 5 The changes in the urban component are shown in Fig. 8 (middle) expressed as NO2 emissions per capita. Recall, that emission rate is the mass divided by the lifetime (that was assumed constant for all areas), therefore, the changes in emissions per capita and the changes in total mass are identical. The relative changes for the urban component (Fig. 8 right)  with the reported emissions: the correlation coefficients between the two data sets from Fig. 9 is 0.9. Note that 2020 US EPA NEI reported emission were incomplete at the time of this study.

The COVID-19 lockdown impact: Europe
The described technique was applied to the European Union countries (plus non-members from former Yugoslavia) where detailed information about the industrial emissions sources is available. The analysis was also done for 3° by 4° areas around 5 36 largest European cities with population greater than 1 million plus some national capitals with population more than 500,000. At the time of this study, the 2020 emissions data were not available, and the 2018-2019 emissions inventory was incomplete with data from many countries missing. Unlike the U.S., the available data were only annual, not monthly, and emissions were typically estimates and not direct measurements. We used the 2017 inventory to obtain the coordinates of the largest sources and then use them in the fitting algorithm. Note, that to avoid double-counting, if more than one city located 10 within an area, we used that area just once (e.g., Manchester and Birmingham are in one area).
The absolute and relative changes between 2018-2019 and 2020 for the three components are shown in Fig. 10. The NO2 decline was particularly large, more than 50%, for the countries in the most western part of the continent where the strictest lockdown measures were taken: France, Spain, and UK. In contrast, the decline in the German, Czech and some other East European cities was only 20-25%. For this reason, two sub-regions were formed for the analysis: Europe-1 (Italy, France, 15 Spain, Portugal, Belgium, Ireland, and UK) and Europe-2 with all other countries. In general, the mean background values and estimated NO2 emissions rates per capita in Europe are similar to those in the U.S. and Canada. However, relative changes are somewhat different.
In 2018-2019, the estimated emissions per capita for both European regions were very similar to those for the U.S. and Canada. In 2020, the urban component declined in almost every analysed area. The average declines for Europe-1 and 20 Europe-2 regions were -52% ±5% and -14% ±8% respectively. This is in general agreement with NOx emission reduction for these two European sub-regions (Guevara et al., 2021). The decline in Europe-1 was rather uniform with all but one area demonstrating a decline of more than 40%. In contrast, only two areas demonstrated a 40% decline in Europe-2, while most of the areas had a decline under 20%. Two areas in Europe-2 (Budapest and Belgrade) demonstrated an increase in NO2. They are located 320 km apart and it is possible that relatively high NO2 values there were cause by some specific meteorological 25 conditions in spring of 2020: the NASA GEOS Composition Forecasting (GEOS-CF) simulations show a positive NO2 anomaly over Hungary in April-May 2020.
As in the case of the U.S. and Canada, the mean background component in Europe shows a smaller decline than the urban component. On average, it was -5%±2% and -12±3% lower in 2020 than in 2018-2019 for Europe-1 and Europe-2 regions respectively, but it was pretty consistent as almost all individual areas demonstrated a decline. Large decline in 30 population-related emissions and relatively small decline in the background component for Europe-1 and the opposite for Europe-2 may create an impression that here is anticorrelation between the background level and population-related component, but it is not true. The large decline in average background for Europe-2 was caused by large negative background https://doi.org/10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. values for the Scandinavian countries in 2020 that also had large negative changes in the urban components. As discussed later in Section 4.5, there is no correlation between the background levels and urban component.
The emissions from industrial sources also demonstrated a decline, although the scattering of the values is large as the changes varies from country to country and from sector to sector. The average decline value is -34%±10 and -13%±16% for Europe-1 and Europe-2 regions, respectively. 5 For illustration purpose, four areas are examined in greater detail below. The Manchester-Birmingham plot (Fig. 11) illustrates a large area of high population density with several power plants to the East. The TROPOMI data analysis shows a 40% decline in the population density-related component and about 25% decline from total emission from the power plants.

The global COVID-19 lockdown impact
To evaluate the COVID-19 lockdown impact worldwide, the analysis described earlier in Section 4 was performed for 263 urban areas around the world. It should be noted that the NO2 "footprints" of cities with the same population vary greatly from region to region with the highest values in Northern Eurasia and Australia, and the smallest in Africa and India.
To illustrate these large differences, Fig. 12 shows the examples of NO2 distribution near cities with population of about 5-6 5 million with very large (Saint Petersburg, Russia) and very small (Dar es Salaam, Tanzania) per capita emissions. The total mass of NO2 per capita related to the urban component for Saint Petersburg was 40 times larger than for Dar es Salaam (all for 2018-2019). While all cities with population greater than 1 million were considered in this study, some of them do not even produce significant NO2 emissions that can be measured by TROPOMI over the three-month period selected for this study. Another obstacle is in Western Africa, where biomass burning made it difficult to estimate "background" levels as they 10 were very different from year to year. In case of China, there are too many cities with population over one million. We raised the limit and considered only cities with population greater than 6 million to keep the number of analysed areas similar to other regions.
The analysis algorithm requires coordinates of individual industrial sources in order to separate them from the urban component. Detailed information about industrial source locations is often not available in many regions. The world power 15 plant database (see section 2.3) was used to locate most of the power plants, while other sources were identified from hotspots on the NO2 maps. A total of 357 such additional sources were identified. Most of them were cement and steel factories, and oil refineries. The complete list of these additional sources is included in the Supplement (file Additional_Sources_World.xls).
In addition, the world busiest airports were included as "industrial" emission sources. However, some sources, in particular, ship tracks may still be missing that may affect estimates for some areas. An example is given in Supplement D. Also, the 20 main highways are not included in the present statistical model. Some of them are identifiable in the residual maps. Ship tracks and highways could be added to the statistical model in the future. Fig. 13 (top). The analysed period from mid-March to mid-June is close to spring in the Northern Hemisphere and autumn in the Southern Hemisphere, i.e. the seasons with very similar values of lifetime (Lange et al., 2021). Therefore, seasonal differences between the two hemispheres should 25 be minimal, and maps of the main estimated components should well represent their global distribution. The highest values are seen over East China and the northern part of Central Europe, while the lowest are mostly over South America and East Africa. In general, it is similar to the global NO2 distribution (Krotkov et al., 2016), although the background levels are much lower than those seen over NO2 hotspots on the mean satellite NO2 map.

The map of mean background component values for all 263 sites is shown in
The values of the urban component, expressed as annual emission rate per capita, are shown in Fig. 13 (bottom). The 30 highest values are over Russia and are likely related to additional NOx emissions due to heating in a relatively cold period in March-April. Another hotspot is Edmonton, but as mentioned, its high value is due to poor separation urban emissions from individual areas for that region and areas with components below some threshold levels were excluded. Fig. 15 shows the maps of changes for individual areas for the background and urban components.
China shows the smallest and not significant decline in the urban component, while the majority of the regions demonstrated statistically significant urban emissions decline within the range 18-28%. The decline was the largest, 36-52%, in three regions: Europe-1, South America and India. The map of the urban emission changes (Fig.15, bottom) shows that the 20 first two regions indeed contain countries with large decline of urban emissions. In case of India, a similar decline can be seen in neighbouring Pakistan and Bangladesh. In Africa, a decline is seen at the south and the north of the continent, while countries in West Africa mostly show no decline and even some increase probably due to a contribution from forest fires.
Changes in the background component are shown in Fig. 14 and Fig. 15 (top). Relative changes of the background component are typically within ±10% and much smaller than in the urban component. One of the exceptions is the 25 Scandinavian countries and Germany that caused an overall -13% change in the background value for the Europe-2 region.
The urban and background components are fairly independent. Analysis of all 263 areas revealed that there is no correlation between changes in urban NO2 emissions and changes in background values (the correlation coefficient is -0.009).
As mentioned in section 4.3, the industrial NO2 component varies from area to area and from one type of NO2 source to another, although there are some clear regional differences. As the main COVID-19 lockdown in China occurred earlier (in 30 February), during the analysed period, Chinese cities demonstrated the smallest changes in both urban and industrial components (-2.8% and +5% respectively). There is, however, one exception. Emissions from Wuhan, the city where the pandemic begun, declined by more than 60%. Industrial emissions also declined, but only by 30%. The background component shows no change there. A very strict Wuhan lockdown ended on April 8, 2020, but during that lockdown, NO2 emissions in https://doi.org/10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License.
Wuhan declined 82% relative to the 2019 level (Ghahremanloo et al., 2021). That strict lockdown period lasted for less than one third of the analysed period, but apparently, it took some time for NO2 emissions to return to the pre-lockdown levels.
The largest industrial emissions decline was observed over the same regions where the largest urban emissions decline was observed: Europe-1 and India. It is likely the severe restrictions during the COVID-19 lockdown period there affected the industrial activity. However, on a larger scale this link is not that obvious. Although the lockdown had impact on industrial 5 sources, the correlation coefficient between changes urban and industrial emissions among all analysed areas is only 0.18.
The uncertainties values in Fig. 14 Table S1). The background component has the smallest variability among the three components typically between 5% and 9%. The urban component variability is between 7% and 17% and the decline observed 10 in the urban component for South America, Europe-1 and India is outside 3-σ limits even for individual areas in these regions.
The industrial component was largely added to separate emission from large industrial sources in the urban areas from urban emissions themselves. Emission from such industrial sources are typically similar or smaller than urban emissions and the variability of the industrial component (10%-30%) is similar or larger than that for the urban component.
To demonstrate that the observed NO2 changes in urban emissions are indeed linked to the restricting measures taken 15 by different countries, the estimated percent NO2 changes in emissions per capita were compared to the Google Each Community Mobility Report data (available from https://www.google.com/covid19/mobility/, accessed on March 1, 2021).
The mobility data represent the changes in the number of people at locations of various type and can be used as a proxy for the urban traffic. The three components of NO2 distribution were compared various characteristics of mobility data. The mean values of the three components of the NO2 distribution (background, urban, and industrial) were calculated for every country 20 by averaging the corresponding values for individual cities. Only countries with two or more cities were used in the comparison. Similarly, the mean characteristics of mobility were calculated for all individual countries. Note that the mobility data were averages of all regions for the entire country, while the NO2 changes were estimated for areas around large cities only. Mobility data for China, Korea and some other countries were not available.
The scatter plot of the mobility and the lockdown-related NO2 changes (Fig. 16) demonstrates a very different 25 relationship between mobility changes and changes in urban and background components. Changes in mobility and urban components are correlated (Fig. 16 left). As expected, the relative changes in the urban component are smaller than the mobility changes as the urban component includes more than just mobility-related traffic. The highest correlation is observed when changes in the NO2 urban component are compared with mobility for "retail and recreation", covering visits to restaurants, cafes, shopping centers, theme parks, museums, libraries, movie theatres, and similar locations. The correlation coefficient 30 between the percent changes in per capita emissions and "retail and recreation" mobility is 0.66 (the probability that there is no correlation is less than 0.0001). The correlations of the urban component with other mobility characteristics are lower (about 0.55 and 0.58 for the "workplaces" and "transit stations" categories, respectively). In contrast, there is no statistically significant correlation (the correlation coefficient is 0.03) between the background NO2 and mobility data (Fig. 16 right). The https://doi.org/10.5194/acp-2021-536 Preprint. Discussion started: 6 July 2021 c Author(s) 2021. CC BY 4.0 License. industrial component shows some correlation (the correlation coefficient is 0.5) with the mobility data, although not as high as for the urban component.
To compare our results with other similar estimates of NOx emission, we looked at a recent study by Lange et al., 2021. In their study, TROPOMI NO2 data were used to estimate emissions from 45 sources worldwide and compared them with the available emissions inventories and some other satellite-based emission estimates. As our study also provides total 5 emissions estimates (the sum of urban and industrial emissions) for 33 of them, we compare our results with those from Lange et al., 2021. While the emission estimation algorithms and approaches are different in the two studies, the results show a good agreement and the correlation coefficient between the two sets of emission estimates is 0.78. As expected, the emission estimated in this study are higher than from Lange et al., 2021. This is because we used larger areas and, typically, there was more than one emission source in the analysed 3° by 4° areas of this study. The details of this comparison are available from 10 Supplement C.
In this study, the population density was used as a proxy for urban emissions. We found that this proxy works well unless the population density was very high and the link between population density and emissions might not be linear anymore. Such exceptions were rare (e.g., some regions of New York City). To test the impact of such nonlinearity, we used a different proxy, the night light data set and all calculations were repeated with them. Unlike the population density, the used 15 night light data reached saturation at highly populated areas and did not grow any further. It was found that the results were similar and the major findings about different reaction of the background and urban components to the COVID-19 lockdown were still valid, although with some differences in absolute values (see section Supplement E). The correlation coefficient between mobility and urban emissions estimated from night light is lower (0.57) that for those from population density suggesting that the population density-based analysis captures the NO2 emissions reduction due to lower traffic better. 20

Discussion and conclusion
Statistical regression analysis was used to separate contribution from industrial sources, urban areas, and background levels to the observed tropospheric NO2 VCD and to study the impact of the COVID-19 lockdown on each component separately. The analysis was done for 263 major urban areas around the world grouped into 13 large regions. The algorithm essentially estimates the total NO2 mass for the three different components that then can be converted to emissions assuming 25 a constant NO2 lifetime (or, more accurately, decay time). A constant value of 3.3 hours was used as the lifetime.
Unlike other similar studies that simply removed the background offset (e.g., (Beirle et al., 2011;Lange et al., 2021)), this study included the background component in the analysis. On a scale of several hundred km (as we analyzed 3° by 4° areas), most of the NO2 mass is typically related to the background component. Even in the areas such as New York City, the background component accounts for 2/3 of the total mass. This explains why the estimated impact of the COVID-19 lockdown 30 in urban areas depends on the size of the analysed area: the bigger the size the more background NO2 it includes and, therefore, the smaller is the NO2 difference between the COVID-19 lockdown and reference periods. Abrupt changes and urban and industrial emissions due to COVID-19 lockdown did not immediately result in a similar decline in the background component. This may explain why large changes in NO2 emissions in urban areas produced a relatively small, about 9% decline in global NO2 (Bray et al., 2021). The importance of background NO2 VCD was previously 10 noted by Qu et al., (2021) and Silvern et al., (2019) when they found that the observed satellite trends in remote areas do not match the expected changes. The background NO2 is anticorrelated with the surface elevation. There are several possible factors that contribute to the background component. It could be related to NO2 in the free troposphere where it may have a much longer lifetime and travel long distances. Satellite measurements are also more sensitive to NO2 in the free troposphere than in the boundary layer and a relatively small amounts of NO2 there produce a larger signal in satellite data. Another possible 15 explanation is that at low concentration in the boundary layer, NO2 may have longer lifetimes than the value of several hours in the plumes. The fact that NO2 fluctuations remain persistent over longer time in clean conditions than over polluted conditions (Vinnikov et al., 2017) indirectly confirms that.
The origins of background NO2 are still largely related to urban and industrial sources as it is clearly higher in the northern hemisphere, particularly over China, Central Europe and Eastern U.S., than in the southern hemisphere and tropics. 20 However, the analysed three-month period may simply be not long enough for the lockdown to cause large changes in the background levels. There are also other NOx sources such as soil emissions (Hudman et al., 2012;Sha et al., 2021). They as well as sources aloft, such as lightning, and to a lesser extent, aircraft NOx directly contribute to the background component.
It is estimated that lightning is responsible for roughly 16% of global production and most of this NOx is found in the free troposphere (Bucsela et al., 2019). 25 The urban and industrial components are based on plume dispersion functions and correspond to NO2 near the ground, probably in the boundary layer. The urban component is based on the population density and the assumption that emissions per capita are uniform everywhere in the analysed area. There are very large differences, up to factor of 40, in estimated emissions per capita between the areas. The estimates were done for 3-month periods. For such short time interval, most of the cities with population more than 1 million produce a statistically significant signal, that can be detected in TROPOMI NO2 30 data. As estimated emissions per capita are rather uniform, they can be used to account for urban component outside large cities. Thus, it should be possible to estimate background, urban, and industrial components on the global scale and analyse the residuals in search of other factors contributing to the NO2 budget.
The approach described in this study can be used to estimate emissions from cities and industrial point sources. For the latter, only source coordinates are required. A comparison of reported and derived from TROPOMI data NO2 emissions for the U.S. demonstrated a good correlation (0.9) between them. As source coordinates can be also detected from satellite data alone (Beirle et al., 2019;Ding et al., 2020;McLinden et al., 2016), it may be possible to develop an independent "topdown" NO2 emission inventory from satellite measurements to complement and emissions improve available "bottom-up" 5 inventories as it was done for SO2 (Liu et al., 2018). This could be important for regions, where no other emission information is available.
This appendix contains additional details of the used fitting algorithm that is largely based on the algorithm for multiple point source emission estimates (Fioletov et al., 2017). TROPOMI NO2 VCD can be expressed as a sum of contributions αi·Ωi from all individual industrial sources (i), a population density-related term αp Ωp, an elevation-related background, and noise (ε): TROPOMI NO2 (θ, φ) = α0 +αp Ωp (θ, φ)+ Σ αiΩi (θ, φ)+ (β0 + β1(θ -θ0) + β2(φ -φ0))·exp(-H(θ, φ)/H0) + ε(θ, φ) (A1) 5 All Ω function are normalized (i.e., their total integral equals to 1) plume functions: the value of that function for a particular pixel with latitude θ and longitude φ, is proportional to the value of the plume parameterization from the source i located at the latitude θi and longitude (φi) (all in radian). The parameterization assumes that the plume is moving downwind along a straight line has a Gaussian shape spread across that line. To describe the plume, we can rotate satellite pixels for a particular day around the source, so the plume would always be moving from north to south, apply the plume parameterization, and then 10 rotate the pixels back. If (xi, yi) and (x′i, y′i ) are the pixel's Cartesian coordinates (km) in the system with the origin at the source i before and after the rotation respectively, then they can be calculated from the pixel and source latitudes and longitudes as: xi= r·(φ-φi)·cos(θi); yi= r·(θ-θi); 15 x′i = xi · cos(-ω) + yi · sin(-ω); y′i = -xi · sin(-ω) + yi · cos(-ω); where r=111.3 km·180/ π (or r=6371 km· π/180 for latitude and longitude in degrees); ω is the pixel wind direction (0 for north); φi and θi are the source i longitude and latitude (all in radian). Note that there was a typo in this original formula for r 20 in Fioletov et al., (2017).
, therefore the parameter αi represents the total observed number of NO2 molecules (or the NO2 mass) near the source i. If TROPOMI NO2 is in DU, and σ is in km, then a is in 2.69·10 26 molec. or 0.021 T(NO2). Furthermore, the emission strength (E) can be calculated as E= α/τ assuming a simple mass balance.
As mentioned in section 3, some of the sources used in the analysis are not point sources but clusters. In that case, Ωi