High-resolution mapping of urban NO₂ concentrations using Retina v2: a case study on data assimilation of surface and satellite observations in Madrid

Urban air pollution poses a significant health risk, with over half the global population living in cities where air quality often exceeds World Health Organization (WHO) guidelines. A comprehensive understanding of local pollution levels is essential for addressing this issue. Recent advancements in low-cost sensors and satellite instruments offer cost-efficient complements to reference stations but integrating these diverse data sources in useful monitoring tools is not straightforward. This study presents the updated Retina v2 algorithm, which generates high-resolution urban air pollution maps by assimilating heterogeneous measurements into a portable urban dispersion model. Tested for NO₂ concentrations in Madrid during March 2019, it shows improved speed and accuracy over its predecessor, with the ability to incorporate satellite data. Retina v2 balances performance with modest computational demands, delivering similar or better results compared to complex dispersion models and machine learning approaches requiring extensive datasets. Using only TROPOMI satellite data, citywide NO₂ simulations show an RMSE of 19.3 µg m⁻³, with better results when hourly in-situ measurements were included. Relying on data of a single ground station can introduce biases, which can be mitigated by incorporating satellite data or multiple ground stations. Including more stations improves accuracy, with 24 stations yielding a correlation of 0.90 and an RMSE of 13.0 µg m⁻³. The benefit of TROPOMI diminishes when data from five or more ground stations is available, but it remains valuable for many cities which have limited monitoring networks.

1 Introduction

More than half of the world's population lives in cities, where most people breath air that exceeds the World Health Organization's (WHO) air quality guidelines (WHO, 2021). Elevated levels of nitrogen dioxide (NO₂), primarily from urban traffic and residential emissions, significantly contribute to this health issue. NO₂ is linked to respiratory diseases, particularly asthma, leading to respiratory symptoms (such as coughing or difficulty breathing), hospital admissions and visits to emergency rooms. According to the WHO air quality database (WHO, 2023), 77 % of the population in the 4000 assessed towns and cities are exposed to mean annual NO₂ levels above the recommended limit of 10 µg m⁻³. This figure rises to 90 % in cities in low- and middle-income countries, paralleling large-scale urbanisation and economic development.

Addressing urban air pollution requires a detailed understanding of local pollution levels. This is best achieved with a dense network of reference stations, as traffic patterns and urban design can cause strong gradients in air pollution (Cummings et al., 2022). However, the high costs of installing and maintaining such networks often leave cities, especially in low- and middle-income countries, without adequate monitoring infrastructure.

In recent years, alternative air quality measurements from low-cost sensors and satellite instruments have become available. The monitoring of urban air quality will greatly benefit from incorporating these complementary measurements (WMO, 2024). However, integrating different data sources in a transparent manner is challenging because they differ in sampling frequencies and spatial representativeness. While low-cost air quality sensors can provide detailed spatial observations in urban areas, they often come with significant uncertainties (Snyder et al., 2013). As satellites observe air pollution of the entire troposphere, the relationship between column concentrations and surface-level concentrations must be resolved first. Also, polar-orbiting satellites pass over the same area only once per day, missing a substantial part of the diurnal cycle (Boersma et al., 2009).

Air quality models are therefore essential to create maps from measurements, as they not only fill in the unsampled areas and times, but also (in the more advanced data fusion methods) consider the different spatial representativity and accuracy of the measurements, and – for satellite measurements – the height distribution in column measurements.

Modelling at the urban scale can be done by Land Use Regression (LUR) models (Hoek et al., 2008), which solve statistical relations between surface concentrations and geographic data. They are commonly used in exposure studies, providing maps at high spatial resolution but lacking a time component. Another approach is to downscale the output of regional chemical transport models to high-resolution sub-grid concentrations (e.g. Denby et al., 2020; Kim et al., 2018). Like LUR models, Gaussian plume models (e.g. listed in Kakosimos et al., 2010) are widely used in urban settings due to their low computational demands. Based on an analytical solution to pollutant transport equations and a detailed emission inventory, they can calculate hourly concentrations of air pollutants at street-level under given meteorological conditions.

Better simulation results are obtained when in-situ measurements are spatially assimilated in modelled concentration fields using kriging techniques (Schneider et al., 2017; Criado et al., 2023) or optimal interpolation (Tilloy et al., 2013; Mijling, 2020). This significantly reduces both local biases and background biases.

The TROPOspheric Monitoring Instrument (TROPOMI) is a nadir-viewing imaging spectrometer aboard ESA's Sentinel-5P satellite. Since May 2018, TROPOMI has provided global observations of air quality from space with an unprecedented spatial resolution (5.6 × 3.6 km² at nadir view since 6 August 2019). This resolution offers coarse information on spatial patterns of air pollution within urban environments. For instance, tropospheric NO₂ column measurements of TROPOMI have been used to estimate NO_x emissions in Paris (Lorente et al., 2019), to predict daily surface NO₂ concentrations in Mexico City (He et al., 2023), and to detect spatiotemporal variations of NO₂ in Madrid (Morillas et al., 2024).

TROPOMI observations are used by Kim et al. (2021) to create hourly NO₂ maps for Switzerland and northern Italy on a 100 m resolution, in a combination with reference measurements, geographical and meteorological data. Fu et al. (2023) also add low-cost sensor data for hourly mapping in Tangshan, China. Both studies show that satellite data can contribute significantly to surface NO₂ mapping, despite its coarse resolution and the fact that it is only available once a day under near cloud-free conditions.

Recurrent complication in urban air quality modelling is the need of an up-to-date emission inventory at a high resolution, and realistic estimations of the regional background concentrations. Also, many urban data fusion applications depend on machine learning (e.g. Kim et al., 2021; He et al., 2023) or a detailed local air quality model (e.g. Schneider et al., 2017; Criado et al., 2023). This complicates portability to other cities where the required input data might not be available.

The Retina algorithm (Mijling, 2020) provides a physics-based and portable approach. It has been developed specifically for observation-based high-resolution modelling of urban air pollution using heterogeneous air quality measurements (i.e. of different accuracy and origin). Central in Retina is the open-source AERMOD dispersion model (Cimorelli et al., 2004), developed by the American Meteorological Society (AMS) and United States Environmental Protection Agency (EPA). The model is driven by meteorology and local emissions constructed from proxy data. Observations are used for emission optimization and spatial concentration assimilation. It generates hourly maps of air pollutant concentrations at street-level.

Retina v2, described in this paper, has undergone significant updates to enhance its speed and accuracy. It is faster and uses less computational resources by using AERMOD only for dispersion kernel calculations (Sect. 2.3.1). The estimation of background concentrations has been improved (Sect. 2.2.1). The NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratios are estimated more accurately by replacing the Ozone Limiting Method with a non-linear regression method (Sect. 2.3.2). The estimation of the sectoral emission factors is better stabilised by implementation of a Kalman filter (Sect. 2.3.4). The spatial assimilation of concentration measurements is improved by including time-dependent dispersion characteristics in the model error covariances (Sect. 2.3.5). Most notably, for the CitySatAir project (part of ESA's EO Science for Society program) we extended the algorithm with an additional module to incorporate tropospheric column concentrations of NO₂ measured with TROPOMI (Sect. 2.3.3).

2 Method

The added value of NO₂ column measurements from space is evaluated through their application in Madrid, Spain, for the period of March 2019. The city's extensive network of NO₂ reference stations allows for the exploration of different measurement configurations, including scenarios with and without TROPOMI observations.

The municipality of Madrid extends over an area of about 40 × 43 km², with a population of approximately 3.4 million people. Urban NO₂ pollution levels are amongst the highest in Europe, regularly exceeding the air quality guidelines set by the WHO (WHO, 2021), which recommend limits of 10 µg m⁻³ for annual averages and 25 µg m⁻³ for daily averages. The city of Madrid ranks first in Europe for mortality linked to NO₂ pollution, according to a recent health impact study by ISGlobal, which analysed nearly 1000 European cities (Khomenko et al., 2022). Traffic and residential emissions are mainly responsible for the high surface concentrations found in the city, as Madrid has no heavy industry or other important NO₂ hot spots in its immediate vicinity.

2.1 NO₂ observations

2.1.1 Reference network

Common equipment to perform reference measurements of NO₂ include the Teledyne API 200E and a Thermo Electron 42i $\mathrm{NO}/{\mathrm{NO}}_x$ analyser, both based on chemiluminescence. A catalytic–reactive converter converts NO₂ in the sample gas to NO, which, along with the NO present in the sample, is reported as NO_x. NO₂ is then calculated as the difference between NO_x and NO. We use an accuracy of 4 % of the NO₂ measurements (GGD, 2023). Note that this might be an underestimation for locations downwind of source areas, as the molybdenum converter also reduces other reactive nitrogen species such as PAN and HNO₃ (especially found in aged plumes) to NO, introducing a positive bias in the NO₂ measurements (Steinbacher et al., 2007).

Madrid has an extensive network of 24 reference sites measuring hourly NO₂ concentrations (see Fig. 1), from which 9 qualify as street stations, 12 as urban background stations, and 3 as suburban background stations. Hourly measurements of the network are published in near real-time as open data at the Madrid Open Data Portal (https://datos.madrid.es, last access: 23 September 2025). Lower concentrations of NO₂ are found in summertime, due to favourable atmospheric conditions (e.g. higher boundary layer height) and less emissions during the holiday period. Highest concentrations are found in wintertime, when monthly averaged values are well above 40 µg m⁻³ at roadside and urban background sites. The network-wide average in March 2019 (36.2 µg m⁻³) closely reflects the annual average (34.5 µg m⁻³).

Figure 1 Reference network for NO₂ measurements in Madrid. The red line indicates the municipality border. Basemap source: © OpenStreetMap contributors © Carto, 2024. Distributed under a Creative Commons BY-SA License.

2.1.2 TROPOMI retrievals

TROPOMI is a nadir-viewing imaging spectrometer aboard ESA's Sentinel-5P satellite. Since May 2018, TROPOMI has provided global observations of air quality from space with an unprecedented spatial resolution (5.6 × 3.6 km² at nadir view since 6 August 2019), enabling it to offer coarse information on spatial patterns and gradients of air pollution within urban environments.

Being in a sun-synchronous orbit at 824 km altitude, S5P overpasses Madrid daily around 13:00 UTC. At times the urban area is sampled from two adjacent orbits, typically around 12:30 and 14:10 UTC. At every overpass the retrieval footprints are located differently, sampling different parts of the urban area.

From the radiance and irradiance spectra the NO₂ slant column density can be derived, which is divided into a stratospheric and tropospheric part. By applying an appropriate air mass factor the tropospheric slant column is converted to a tropospheric vertical column density and its accompanying averaging kernel. We use the latest reprocessed product for TROPOMI tropospheric NO₂ columns, version 2.4 (Eskes et al., 2022; van Geffen et al., 2022), which implements a new surface albedo climatology derived from TROPOMI observations, including the viewing-angle dependence of the scattering at the surface.

Although it is recommended that for straightforward application only retrievals with a quality value ≥ 0.75 should be used (i.e. valid retrievals with cloud radiance fractions below 50 %), this criterion is relaxed to ≥ 0.5 (i.e. valid retrievals, including under cloudy conditions) as the averaging kernel is carefully applied (see van Geffen et al., 2022). Also, only footprints which cover the studied domain by more than 50 % are used. In this way, there are on average 14.2 valid retrievals found per day in March 2019.

2.2 Model input data

2.2.1 Background concentrations

An important fraction of the air pollution is transported from upwind regions (see e.g. Harrison, 2018). Using realistic background concentrations of NO₂ is therefore crucial as the dispersion simulation only accounts for locally generated NO₂ within the domain.

Here we use hourly data from the European air quality ensemble from the Copernicus Atmosphere Monitoring Service (CAMS) (Marécal et al., 2015), which for 2019 data consist of 9 state-of-the-art numerical air quality models. More specifically, the ensemble median of the validated reanalysis is used, a data product for which each ensemble member assimilated validated hourly observations of air pollutants as reported to the European Environment Agency. The CAMS regional ensemble covers Europe at a resolution of 0.1 × 0.1°. Note that the coarseness of the data product makes it unsuitable for monitoring urban air quality in detail, as the strong concentration gradients found around strong sources will be averaged out, leading to an underestimation of NO₂ concentrations, see Fig. 2.

Figure 2 NO₂ surface concentrations in Madrid, March 2019, from bilinear interpolation of the average values of the CAMS reanalysis. CAMS grid cells are shown in thin black lines. Also shown are the mean concentrations measured by the reference network.

To avoid double counting of locally produced NO₂ concentrations, we take for each hour a weighted average of the CAMS concentrations found along the municipal perimeter. Let b_C(x) represent the interpolated CAMS concentration at location x, e_v the unit vector along the wind direction, n the normal vector on the perimeter (pointing outwards), and L the part of the perimeter where ${\mathit{e}}_{\mathrm{v}}\cdot \mathit{n}<\mathrm{0}$ (i.e. where background pollution is flowing into the municipal domain). A uniform wind direction is taken, based on the 10 m wind at the domain centre. Then the weighted average for this hour is calculated from the line integral

$$\begin{array}{}\text{(1)}& b={\displaystyle \frac{\mathrm{1}}{W}}{\int }_L{b}_{\mathrm{C}}\left(\mathit{x}\right){\mathit{e}}_{\mathrm{v}}\cdot \mathit{n}dl,\mathrm{with}W={\int }_L{\mathit{e}}_{\mathrm{v}}\cdot \mathit{n}dl,\end{array}$$

where dl represents the elementary arc length along the perimeter. The resulting background value b is taken homogeneously for the entire domain, thus avoiding the use of more advanced advection schemes. A validation of this method for Madrid can be found in Sect. S1 in the Supplement. The background concentration of ozone, needed in the calculation of the NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratio, is calculated from CAMS ozone fields in the same manner.

2.2.2 Meteorology

Meteorology is an important ingredient for dispersion modelling, as it determines how the air pollutant is transported horizontally and vertically. AERMET, the meteorological preprocessor of AERMOD, requires both surface meteorological data (cloud cover, temperature, humidity, dew point temperature, pressure, precipitation, wind speed and direction) and upper air meteorological data (temperature, humidity, wind speed and direction in vertical layers). The wind speed and direction are also used to determine the influx of background concentrations and the main axes of the model error covariance (Sect. 2.2.1 and 2.3.5).

We use the collection of short-range forecasts (issued at 00:00 and 12:00 UTC) from the archive of the European Centre for Medium-Range Weather Forecasts (ECMWF) at the supercomputing facility in Bologna. It is retrieved as 3-hourly output at a 0.05° spatial resolution. Hourly meteorological fields are obtained by temporal interpolation and then interpolated to a representative location central in the Retina domain.

2.2.3 Emission proxies

An important unknown when modelling street-level urban air pollution is the location and strength of the urban emissions. For many cities this information is either unavailable or outdated. By describing traffic and residential emissions with proxies taken from open data sets (see Fig. 3) we enable a versatile model setup which can be applied easily to other cities. The domain boundary is extended outward by 1500 m to account for contributions from sources close to the municipal border. Other sectoral emissions, e.g. from industry, will be accounted for indirectly in either an increased background field or in additional residential emissions.

Figure 3 Proxies used for the estimation of urban emissions. The Retina algorithm distinguishes between highways, primary roads, and residential emissions. Nearby emissions outside the municipal border are also considered.

Road location and road type classification is taken from OpenStreetMap (OSM). All road segments labelled “motorway” or “trunk” are linked to the highway class. All “primary”, “secondary” and “tertiary” segments are linked to the primary road class.

The Madrid Open Data Portal publishes vehicle counts at approximately 4000 locations, mainly from inductive loop sensors at traffic junctions. The data files contain vehicle counts in 15 min bins, which are aggregated into one-hour bins to align with the temporal resolution of the simulations. For each location, a monthly averaged traffic volume cycle is calculated, separated by hour of day and day of week.

Between counting locations, the traffic flow is estimated by spatial interpolation using inverse-distance weighting. The interpolation is done separately for vehicle counts at highways and primary road networks, as they have incomparable traffic volumes.

Population density is a good proxy for residential emissions from activities such as heating and cooking. We use the population densities from the Global Human Settlement project (Freire et al., 2016) which are provided on a 250 m resolution. To reflect the observation that per capita residential emissions decrease when people live closer to each other (e.g. Makido et al., 2012), the emission fluxes are scaled proportionally to the square root of the population density.

The proxy data P_s for sector s (here traffic and residential) are gridded on a high-resolution 10 m grid to enable fast application of the dispersion kernel (see Sect. 2.3.1). Sectoral NO_x emissions for a grid cell indexed with (i,j) are calculated by applying the emission factors f_s to the proxy data:

$$\begin{array}{}\text{(2)}& E_s(i,j)=f_s{P}_s(i,j)\end{array}$$

Emissions change over the day. For traffic emissions this is described by a time dependency in the proxy data. The diurnal cycle in residential emissions, however, is described by 24 different emission factors (each for one hour), as its proxy data are constant in time.

2.3 The revised Retina algorithm

Retina uses past observations for emission optimisation (minimizing the general model bias) and current observations for spatial concentration assimilation (reducing local model biases).

In the emission optimisation step (represented in Fig. 4) the algorithm fits emissions factors for the proxy emissions that best match the spatio-temporal concentration patterns observed by the reference network. The optimisation is repeated every 24 h (see Sect. 2.3.4), using the emission factors and covariances of the previous estimation as a priori. This approach avoids the need of detailed knowledge of vehicle fleet composition and solves mismatches between theoretical and real-world emissions. It also compensates for model biases resulting from incorrect chemistry (e.g. lifetime) and unaccounted (seasonal) emission cycles.

Figure 4 Schematic representation of Retina's emission optimisation workflow. Starting with a priori emission factors, the AERMOD dispersion model simulates NO₂ concentrations at the observation locations for a 24 h period. The Kalman filter infers from the difference between observation and simulation the best update for the emission factors. These values are passed to the next analysis period.

With the most recent estimation of emission factors, Retina simulates the surface concentrations for a specific hour at all 62 266 locations of a non-regular, road-following mesh (see e.g. Lefebvre et al., 2011).

Next, in-situ observations are spatially assimilated in the simulated concentration field using optimal interpolation (Daley, 1991), see Fig. 5. This technique allows for the assimilation of surface measurements with different error margins. At the observation locations, model values are corrected towards the observations. In the surrounding areas, the balance between the model and observation errors determines how the simulation is adjusted (see Sect. 2.3.5).

Figure 5 Schematic representation of Retina's hourly simulation and assimilation workflow. The dispersion model uses optimised emissions and background concentrations to simulate surface concentrations. In-situ observations (if available) are assimilated in an optimal interpolation scheme.

2.3.1 Dispersion kernel

Due to the large number of emission sources, using a straightforward AERMOD configuration can result in long calculation times, especially if simulations must be performed at several vertical levels to recreate the tropospheric columns. Instead, we adopt the approach suggested by Masey et al. (2018), using AERMOD exclusively to calculate the dispersion kernel. This is the dispersion of a unit NO_x emission for a specific hour under given meteorological conditions (e.g., wind speed and direction, atmospheric stability, and boundary layer height).

Removal processes of NO_x are modelled with an exponential decay, which is included in the kernel calculation. In urban settings, the typical lifetime of NO_x is on the order of a few hours (Beirle et al., 2011), and changes significantly when plumes travel from source areas and mix with clean air (Krol et al., 2024). However, as emitted NO_x has a relatively short residence time in urban areas, the specific value of its lifetime is not very critical. We use a heuristic value of 2 h throughout this study.

Dispersion kernels are computed for all emission release heights of the emissions and all receptor heights. Sectoral release heights are 0.5 m for traffic and 10 m for residential emissions. These sector-specific dispersion kernels K_s are then gridded onto the regular high-resolution grid, aligning with the emission grid E. Examples of dispersion kernels are illustrated in Fig. 6.

Figure 6 Examples of dispersion kernels at different release heights and winds (i.e. different atmospheric stability) for surface concentration calculation. The receptor height of the kernels is at 1.5 m. Distances are given in meters.

The NO_x concentration $c^{{\mathrm{NO}}_x}$ (in NO₂ mass equivalent per volume) for a receptor located in grid cell (i ,j) can now be calculated as a superposition of all dispersed emissions, for all contributing emission sectors s

$$\begin{array}{}\text{(3)}& c^{{\mathrm{NO}}_x}\left(i,j\right)={\sum }_s{\sum }_{i^{\prime }{j}^{\prime }}{K}_s(i-i^{\prime },j-j^{\prime })E_s(i^{\prime },j^{\prime })\end{array}$$

This can be interpreted as an element-wise matrix multiplication (i.e., the Hadamard product) of the mirrored dispersion kernel with its origin at (i,j) with the entire emission grid. In first order, the NO_x concentrations depend linearly on contributing sources; an NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratio is applied afterwards (see Sect. 2.3.2). Transport over longer distances must be accounted for to prevent underestimation: although a single grid cell at a long upwind distance may contribute only a small amount, the cumulative effect from numerous such grid cells becomes significant. In our algorithm we consider contributions up to 30 km to account for transport effects in the entire simulation domain. To ensure computational efficiency, contributions at longer distances are computed at lower resolutions, as emission source locations become less critical with increasing distance.

2.3.2 Surface concentration simulation

We calculate NO₂ concentrations from NO_x concentrations using a time and location dependent NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratio. The dynamic equilibrium between NO and NO₂ is affected by temperature, available ozone (which generates NO₂ from NO), and solar radiation (which generates NO from NO₂). The local NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratio is hourly estimated using parameters available during simulation:

• the local simulated NO_x concentration.

• the background O₃ concentration, taken from the regional CAMS ensemble.

• the background NO₂ concentration, also taken from the regional CAMS ensemble.

• the temperature, as a measure of reaction speed for conversion NO to NO₂.

• the Solar Elevation Angle (SEA), as a measure of radiation available for photolysis of NO₂.

As the dependence on these parameters is non-linear, we train an extreme gradient boosting (XGBoost) model (Chen and Guestrin, 2016) to estimate NO₂ ratios at simulation time (see Sect. S2).

From the NO_x concentration $c^{{\mathrm{NO}}_x}$ calculated by Eq. (3), the NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratio r, and the background concentration b (from Eq. 1), we calculate the NO₂ surface concentration c as

$$\begin{array}{}\text{(4)}& c(i,j)=b+r(i,j)c^{{\mathrm{NO}}_x}(i,j)\end{array}$$

2.3.3 Column concentration simulation

For column concentration simulations, which are necessary for comparison against satellite observations, we use the same approach as for surface simulations but with different settings (see Table 1).

Table 1 Overview of simulation settings.

Given the large footprints of the satellite observations relative to the urban domain, it is crucial to maximize the information content from single NO₂ retrievals. Gridding and averaging to a model grid, or clustering orbits in time, would result in valuable information loss. Instead, we project simulation results to individual footprints in individual orbits, following these steps:

1. Modelling of NO_x concentrations at high resolution at 9 different heights. We use the vertical grid definition from the CAMS regional ensemble from surface to 5 km.

2. Spatially aggregating the simulated values to individual satellite footprints. The coarser resolution in the horizontal (250 m) is justified as the satellite footprints are in the order of kilometres.

3. Temporal interpolation of simulation values to the exact satellite overpass time.

4. Applying the averaging kernel associated with the retrieval, after conversion of the concentration profile to partial columns matching the kernel's levels. This reduces errors resulting from profile assumptions in the retrieval method (Eskes and Boersma, 2003).

Finding a realistic NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratio for column simulation is not straightforward, as the chemical equilibrium changes with height. Close to sources, lower layers exhibit smaller ratios at low temperatures (when conversion from primarily emitted NO is not established yet), while higher layers show reduced NO₂ due to stronger photodissociation from solar irradiance. An intricate chemical analysis is avoided by taking the ratio R from columns of the CAMS regional ensemble (see Sect. S3). For Madrid, R fluctuates between 0.4 and 0.8, with the lowest values found in winter. Note that NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratios in columns are generally higher than surface ratios due to increased ozone availability.

The column simulation includes only local contributions, excluding the background column (i.e., concentrations from emissions outside the domain and the free tropospheric column above the boundary layer). To simplify matters, we do not simulate this background column. Instead, the background column concentration is determined for each overpass by fitting the simulated column concentrations to the observations.

Let $C_{kls}^{{\mathrm{NO}}_x}$ be the simulation of the partial NO_x column concentration for the footprint of retrieval k, for atmospheric layer l, and for emitting sector s. $C_{kls}^{{\mathrm{NO}}_x}$ is calculated from dispersed proxy data following step 1 to 3. Let A_kl be the averaging kernel element for layer l and retrieval k. With the ratio R known from CAMS and the background column B fitted, we can write for the simulated tropospheric NO₂ column C:

$$\begin{array}{}\text{(5)}& C_k=B+R_k{\sum }_l{A}_{kl}{\sum }_s{C}_{kls}^{{\mathrm{NO}}_x}\end{array}$$

Figure 7 demonstrates that this approach effectively simulates the location, shape, and strength of the pollution plume over the city. The y-intercept in the scatter plots represents the estimated background concentration B, which is determined using ordinary least squares regression, using the reciprocal of the retrieval error as error weights.

Figure 7 Three examples of tropospheric NO₂ columns over Madrid under different wind speeds and directions (indicated by the black arrows). Panels (a), (d), and (g) show retrievals from TROPOMI in single overpasses. Panels (b), (e), and (h) show corresponding column simulations (including background columns) by Retina, where emission factors have been optimised before against surface concentrations. The background columns are estimated from the linear fit between simulations from local emissions and observations (c, f, i).

2.3.4 Emission optimisation: estimating emission factors

In the dispersion modelling described above, the sector-specific emission factors (f_s) remain unknown. These can be estimated from the observations. Rewriting Eqs. (2) to (5), we can express simulations of both NO₂ surface concentrations at location i and tropospheric NO₂ column concentrations at footprint k as a linear combination of f_s:

$$\begin{array}{}\text{(6)}& c_i=b+r_i{\sum }_s{\mathit{\alpha }}_{is}{f}_s\mathrm{and}{C}_k=B+R_k{\sum }_s{\mathit{\beta }}_{ks}{f}_s\end{array}$$

Here, α_is and β_ks are calculated from the dispersion of sectoral proxy data by the dispersion kernels, b and B are the background concentration at the surface and the background column concentration respectively, and r and R the corresponding NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ ratios.

Finding the emission factors from observations is an ill-posed inverse problem, which we regulate here with a Kalman filter. The technical description of the estimation is described in Sect. S4. For in-situ observations, the state vector consists of 25 unknown f_s values: one emission factor for traffic, and 24 elements describing the diurnal cycle of the residential emissions. This is different when using TROPOMI observations only, as they only provide information around overpass time. In this case, we use an a priori diurnal cycle for residential emissions (see Sect. S5). The state vector is then reduced to two unknowns: one emission factor for traffic and one factor that scales the a priori residential diurnal cycle.

By carefully selecting the covariances in the filter we optimise the response time without introducing too much noise from measurement and model errors. Starting with an arbitrary state vector, a spin-up time of at least one month (i.e. ∼ 30 optimisation iterations) is needed. Figure 8 shows an example how satellite observations of a single orbit over Madrid are used for an updated estimation of the emission factors.

Figure 8 Example of the Kalman filter applied to TROPOMI observations in a single orbit over Madrid on 8 March 2019. The grey dots represent the simulated tropospheric column concentration values, which show a clear underestimation for that day. After applying the Kalman filter the emission factors are updated, resulting in new simulations (blue dots) that are closer to the x=y line. The updated emission factors are used as a priori in the emission optimisation of the following day.

2.3.5 Spatial assimilation of surface concentrations

Retina uses the optimised emissions to simulate hourly surface concentration fields, which serves as a priori for the spatial assimilation of in-situ measurements. This assimilation process corrects local model errors from e.g. incorrect proxy data or inhomogeneities in the background field. With vector x_f representing the simulated surface concentration field (i.e. the forecast), and vector z containing the in-situ observations, the statistical interpolation can be written as:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{x}}_{\mathrm{a}}={\mathit{x}}_{\mathrm{f}}+\mathbf{K}(\mathit{z}-H\left({\mathit{x}}_{\mathrm{f}}\right))\text{(8)}& {\displaystyle }\mathbf{K}={\mathbf{P}}^{\mathrm{f}}{\mathbf{H}}^{\mathrm{T}}\left({\mathbf{HP}}^{\mathrm{f}}{\mathbf{H}}^{\mathrm{T}}+\mathbf{R}\right)\end{array}$$

The update of the forecasted field (x_a−x_f) depends on the difference between the observations z and the collocated simulations H(x_f). Matrix H maps the simulations to the observation locations, and R contains the observation covariances, as in Mijling (2020). The update is determined by the Kalman gain matrix K which balances between the observation error and the model error. Figure 9 illustrates this spatial assimilation cycle.

Figure 9 Spatial assimilation for 7 March 2019, 13:00. Starting from the top-left panel and moving counter-clockwise: (a) in-situ observations and the forecasted field based on optimised emission factors; (c) difference between observations and simulation; (d) concentration update using optimal interpolation; (b) a posteriori field with assimilated observations.

The model error covariance matrix, P^f, represents the spatial extent of model errors: an error at the observation location implies errors in nearby areas. This covariance is approximated by accounting for the spatial representativity of observations, which varies between street and background locations, and by incorporating hourly changes in atmospheric dispersion, writing:

$$\begin{array}{}\text{(9)}& \mathrm{cov}\left({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}},t\right)={\mathit{\sigma }}_{\mathrm{1}}{\mathit{\rho }}_A\left({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}}\right){\mathit{\rho }}_B\left({\mathit{x}}_{\mathrm{1}}-{\mathit{x}}_{\mathrm{2}},t\right){\mathit{\sigma }}_{\mathrm{2}},\end{array}$$

where σ₁ and σ₂ are the model errors at location x₁ and x₂. ρ_A represents the correlation between modelled time series at these locations (which is here calculated from all simulations for March 2019). Correlations between traffic locations and background locations will be lower than between similar locations, therefore ρ_A reduces the impact of an observation on areas it is not representative of (see Fig. S7).

New in Retina v2 is the inclusion of ρ_B, representing the spatial correlation due to the time-dependence in pollutant dispersion. We want to describe this in terms of the two-dimensional dispersion kernel introduced in Sect. 2.3.1. As demonstrated in Sect. S6, ρ_B can be calculated by multiplying the 2D kernel with a copy of itself, shifted by the distance x₁−x₂. This results in a dispersion correlation field which is symmetric along the downwind and crosswind axes. The correlation at a distance d along each main axis can be approximated with

$$\begin{array}{}\text{(10)}& {\mathit{\rho }}_B\left(d\right)\approx {\left(\mathrm{1}+\left|d/L\right|\right)}^{\mathrm{0.75}}\exp \left(-{\left|d/L\right|}^{\mathrm{0.75}}\right),\end{array}$$

with L the fitted correlation length for the considered hour. Typically, the correlation lengths are larger at night and shorter during the day. The exponent 0.75, determined heuristically, provides representative high correlations around the measurement location and best captures the decaying tail.

3 Results

3.1 Model accuracy with TROPOMI-only data

As explained in Sect. 2.3.3 and 2.3.4, retrievals in individual overpasses of TROPOMI can be used to optimise the emissions for the dispersion model. We evaluate the results when optimisation is done at 24 h intervals, with a two-month spin-up time to ensure convergence from a priori values. The space-derived emission factors are used to simulate hourly surface concentrations at a high resolution, as shown in Fig. 10.

Figure 10 Simulation of NO₂ surface concentrations by Retina averaged for March 2019, using TROPOMI observations for daily estimates of emission factors.

The hourly NO₂ concentrations are validated against time series at the 24 sites of the reference network. The statistics of the hourly time series of Retina + TROPOMI and the regional CAMS ensemble are listed in Table 2. We use three different statistical parameters for evaluation:

• Correlation. the Pearson correlation coefficient. A value closer to 1 indicates a better capacity to describe the hourly dynamics of NO₂ concentrations. However, the model might still be off by a bias and/or a multiplication factor.

• Bias. difference between the average simulated value and the average observed value in the considered period. A negative value indicates a systematic underestimation of the simulation, while a positive value indicates a systematic overestimation of the simulation.

• RMSE. the root mean square error. A lower value indicates a smaller distribution of simulations around the true values. Note that the RMSE can be dominated by a bias.

As expected, the interpolated CAMS reanalysis shows a strong negative bias (−10.0 µg m⁻³ on average), particularly at roadside locations. The correlation, in contrast, is quite high (0.86 averaged over all 24 reference sites). This can be attributed to the CAMS ensemble members assimilating observations from a selection of background stations, namely ES0124, ES0126, ES1532, ES1939, ES1942, ES1947, and ES1946 (see CAMS, 2022).

Table 2 Retina validation statistics for hourly surface NO₂ concentrations when optimising emissions using TROPOMI data (validation against interpolated CAMS concentrations in parentheses). Average statistics are indicated in bold font.

Dispersion modelling based on TROPOMI-estimated emissions produces realistic ground concentrations: the absolute biases are reduced at most reference sites (from 10.0 to 0.8 µg m⁻³ on average), particularly at roadside and urban background locations. However, this improvement comes at the cost of introducing more scatter, resulting in lower correlation (0.74 on average), and an RMSE increasing from 18.1 to 19.3 µg m⁻³.

Satellite observations can only estimate emissions around overpass time. Wrong assumptions in the diurnal cycle for other hours introduce an additional error. Therefore, the statistics improve when evaluated for overpass hours only (12:00 to 14:00 UTC): the correlation is higher (0.81) and RMSE is lower (10.9 µg m⁻³).

A large negative bias is still found at ES0118 (Escuelas Aguirre), which is close to a busy intersection. Retina might underestimate the local traffic flow here, or the additional pollution burden due to deceleration and acceleration of congested traffic. A notable overestimation of NO₂ concentrations occurs at location ES1525 (Cuatro Caminos). This is subject to further investigation; it might be related to an overestimation of local traffic intensities.

Section S7 presents the city-wide performance for all months in 2019. March 2019 is approximately representative of the annual average, although seasonal variability is evident. The highest correlation (0.86) is observed during the colder months, when NO₂ concentrations are elevated. In summer, the correlation decreases to 0.70, but the RMSE is lower due to reduced bias and overall lower NO₂ levels.

3.2 Model accuracy under different network configurations

The Retina algorithm can also ingest in-situ measurements of one or more ground stations to estimate the emissions factors. The extent and the distribution of this observation network will affect the quality of the emission estimations and therefore the accuracy of the NO₂ simulations.

First, we assess the influence of station location on emission optimisation using a single reference station, both with and without satellite data. Table 3 summarises the average validation statistics for hourly simulated NO₂ concentrations across all 24 observation locations.

Table 3 Mean validation statistics using different reference stations for emission optimisation, sorted from high to low bias. The right column lists the corresponding total emission estimation. Outside parentheses: including TROPOMI data; inside parenthesis: excluding TROPOMI.

The emission optimisation depends on the selected location. Local model issues or biased reference measurement can lead to biases at other locations. For example, using data of road station ES1943 (Plaza Elíptica) introduces an overall bias of 17.2 µg m⁻³. Retina apparently underestimates the high concentrations at this site, and compensates by adding emissions, resulting in overestimations elsewhere. Combining the data with satellite data reduces the bias significantly. This improvement can be seen at all sites with large absolute bias. The impact of the satellite measurements remains limited, however, as only corrections around overpass time can be provided.

For 13 stations the mean absolute bias is below 4 µg m⁻³. At these sites Retina realistically describes the local concentrations with the provided traffic and residential proxies. Using one of these stations for emission optimisation, the RMSE and correlations are better than when using TROPOMI measurements alone, as the in-situ data provide valuable information on the entire diurnal emission cycle. As can be seen from the table, adding satellite data to these ground data does not improve the results significantly.

Increasing the number of in-situ stations enhances the accuracy of simulation results due to compensating errors, especially when the stations represent a balanced mix of street and background locations. The table shows the results of a test with 5 stations: ES1525, ES1940 (both street stations); ES0126, ES1937 (both background stations); and ES1947 (a suburban background station). It can also be seen that the effect of adding TROPOMI data becomes negligible, as the amount of valid daily satellite retrievals is small compared to the 120 daily measurements made by the 5 stations.

A final test incorporates all 24 reference stations for emission optimisation. No significant change in validation statistics is observed compared to the 5-station scenario. The resulting RMSE, 17.2 µg m⁻³, can thus be regarded as the systematic error inherent in the Retina approach for hourly NO₂ simulation.

Compared to this full in-situ analysis, TROPOMI-based emissions tend to attribute more emissions to traffic and less to residential activity (see Fig. S10), resulting in up to 5 µg m⁻³ higher concentrations on roads and up to 1.5 µg m⁻³ lower concentrations in urban backgrounds.

3.3 Results of spatial concentration assimilation

Spatially assimilating the in-situ observations as described in Sect. 2.3.5 reduces simulation biases in the vicinity of an observation location (e.g. due to wrong local emissions), while at longer distances it reduces simulation errors due to inaccuracies in hourly background concentrations. Table 4 compares the validation results before (i.e. the plain dispersion simulation) and after the spatial concentration assimilation of reference measurements. We use a leave-one-out cross-validation: at each validation location the observations of the other 23 locations are assimilated in the simulation fields. Based on the average NO₂ concentrations (36.2 µg m⁻³) and the average RMSE found in the validation (17.2 µg m⁻³), the relative error for simulation of hourly NO₂ concentrations is estimated to be 48 %. The data assimilation increases the correlation from 0.79 to 0.90 and reduces the RMSE to 13.0 µg m⁻³. Therefore, the relative error improves to 36 % after assimilation, with local systematic biases remaining as the primary source of error.

Table 4 Leave-one-out cross validation statistics for spatial concentration assimilation using all reference stations. (Statistics for model simulation, based on 24-station emission optimisation, inside parentheses). Average statistics are indicated in bold font.

^a Distance to the nearest observation site (in km), ^b Number of measurements, ^c Mean observation at this location, ^d In units of µg m⁻³.

Note that an additional bias can be introduced at places where the covariance is wrongly defined. This happens for instance at street location ES0118 and city park ES1939 (El Retiro), which are at 800 m distance. The high concentrations found at ES0118 influence the spatial interpolation towards the urban background. Vice versa, the low concentrations measured in El Retiro park propagate towards the nearby street site, contributing to a negative bias.

Figure 11 illustrates the performance of Retina at a street location, an urban background location, and a suburban background location. It shows the hourly NO₂ series for a representative week when only satellite observations are used in the simulation, representing the case for a city without any air quality observation network. This is compared with time series where all in-situ observations are used, representing a city with an extensive ground network. Not surprisingly, the best results are obtained when all reference measurements are used.

Figure 11 Comparison of NO₂ time series at different locations for a week in March 2019. The red line represents simulations from the Retina algorithm using only TROPOMI observations for emission optimisation; the blue line represents the leave-one-out time series using data from all other reference stations.

Figure 12 shows the average NO₂ concentration map of Madrid based on all hourly concentration fields of March 2019. Highest concentrations are found on and near the highways, such as the M30 surrounding the city centre in the East and South. Lowest concentrations are found in the sparsely populated El Pardo area in the north. Local concentration reductions are found in e.g. El Retiro park. Note the accumulation of air pollution in the southwest area of the municipality due to predominant winds (1–3 Beaufort) from the northeast in this period.

Figure 12 Average surface NO₂ concentrations in Madrid for March 2019, using the Retina algorithm with hourly data from 24 reference stations for emission estimation and spatial assimilation. The right panel zooms in on an 5 × 5 km² area around Retiro city park, where concentrations are notably lower than in nearby roads and built-up areas.

4 Discussion

Retina v2 introduces a more realistic NO₂ chemistry scheme, an improved background concentration estimation, and a better stabilised emission optimisation. The spatial assimilation of measurements is improved by including time-dependent dispersion characteristics in the model error covariances. Most notably, the algorithm is now capable to incorporate tropospheric column concentrations retrieved with satellite instruments, such as TROPOMI on the Sentinel-5 Precursor (S5P) satellite.

Satellites in polar orbits, such as S5P, pass over the same area just once each day, therefore missing a substantial part of the diurnal cycle of NO₂. As a result, directly assimilating NO₂ satellite observations into concentration fields has limited utility, given the short atmospheric lifetime of NO₂ which limits the system's memory to just a few hours. Instead, we use the satellite retrievals to improve estimations of urban NO_x emissions. As the number of daily TROPOMI observations over the urban area is limited (14 on average for the Madrid domain), it is important to get the most out of each satellite retrieval by interpolating the model simulations to the individual footprints at exact overpass time. Applying the averaging kernel minimizes errors resulting from profile assumptions in the retrieval method.

4.1 Comparison with relevant studies

NO₂ concentrations in urban areas vary strongly in space and time. Unsurprisingly, the CAMS regional ensemble is unfit to represent local NO₂ concentrations in urban areas. Due to its coarse resolution, its interpolated values underestimates concentrations by 10.0 µg m⁻³ in Madrid in March 2019. However, the CAMS ensemble provides valuable input data for background concentration estimation and NO${}_{\mathrm{2}}/{\mathrm{NO}}_x$ column ratios for downscaling algorithms such as Retina.

The validation results show that an urban dispersion model can successfully be built based on CAMS input data and proxy data for traffic and residential emissions. Validation of the hourly NO₂ simulations based on periodic emission optimisation by TROPOMI show a reduction of the mean bias to 0.8 µg m⁻³ and an average RMSE of 19.3 µg m⁻³. Part of the error is caused by wrong assumptions in the diurnal emission cycle, as TROPOMI is only able to capture the emissions around its overpass time. Finding a better a priori diurnal and weekly emission cycle is subject to further investigation.

CALIOPE-Urban, an advanced dispersion model based on a detailed emission inventory and running on a supercomputer, produces an RMSE of 23 µg m⁻³ for hourly simulations in 2019 for the city of Barcelona (Criado et al., 2023). This RMSE reduces to 16 µg m⁻³ when reference data of 12 stations is spatially assimilated using Universal Kriging, and further to 12 µg m⁻³ when also an additional basemap layer based on Palmes-tube measurements is included. Schneider et al. (2017) find a citywide RMSE of 14.3 µg m⁻³ for a similar data fusion method of 24 low-cost sensors in the EPISODE model for Oslo in January 2016. From Table 4 can be seen that these figures are comparable to the RMSE of 13 µg m⁻³ by Retina when all reference measurements are spatially assimilated.

Alternatively, several studies use a machine learning approach to generate hourly surface concentrations maps from a collection of data sets. While our study focuses specifically on the urban area, these approaches typically cover larger regions and incorporate a broader and more diverse set of in-situ measurement locations. Kim et al. (2021) train a predictive model including data from TROPOMI and 340 reference stations in Switzerland and northern Italy, resulting in a spatio-temporal correlation of 0.79 with 40 test sites. Tables 3 and 4 show that the correlation of Retina simulations with reference measurements in Madrid is 0.74 when only TROPOMI is used (i.e. no surface measurements), increasing to 0.79 when 5 or more surface stations are also used for emission optimisation. Best correlation (0.90) is obtained when the reference measurements are also spatially assimilated.

Fu et al. (2023) use data from 266 reference stations and 666 low-cost sensors (LCS) in the Tangshan area (East China). TROPOMI data are used in XGBoost models to fill in missing data at reference sites and to enhance the observations at LCS sites. When trained with reference data only, their predictive model has a correlation of 0.79 and an RMSE of 17.1 µg m⁻³ . The RMSE improves to 16.9 µg m⁻³ when including TROPOMI, and further to 16.3 µg m⁻³ when all LCS observations are also included in the training. Table 3 shows that for Retina-Madrid the RMSE is 17.0 µg m⁻³ when 5 reference stations for emission optimisation are used.

By adding more in-situ data, the RMSE of hourly simulations remains around 17 µg m⁻³, corresponding to a relative error of 48 %, which can be considered as the systematic error of the Retina dispersion modelling. This is an improvement over the previous version described in Mijling (2020), which had an estimated error of 58 %. More research is needed to further reduce this error by addressing its various sources, such as:

• Using better proxy data, particularly a more realistic traffic model for the relative distribution of traffic volumes.

• Including emission hot spots from industry and power generation.

• Improving local dispersion modelling by accounting for e.g. traffic junctions, speed bumps, and street canyons.

• Using more realistic estimates of the background concentration field.

• Improving NO_x chemistry, e.g. by introducing variable NO₂ lifetimes.

Urban NO_x emissions can be calculated from the emission optimisation results by summing Eq. (2) over hours and grid cells. The observation-based NO_x emission estimates for March 2019 in Madrid vary between 1185 Mg NO (when derived from TROPOMI observations) and 1336 Mg NO (when derived from in-situ observations at 24 locations). This is larger than the 705 Mg NO found in the CAMS emission inventory for this month (Soulie et al., 2024) but is in correspondence with the DECSO v6.3 inventory (van der A et al., 2024; based on TROPOMI observations and chemical transport model calculations) being a factor 2 larger than the CAMS inventory. The potential of Retina to estimate realistic urban emissions is subject of further investigation.

4.2 Calculation time

The Retina v2 algorithm is implemented in Python scripts. Calculations were performed on a Linux workstation with an Intel Core i7-9700 at 3GHz, having a single CPU and 8 cores. The total calculation for a high-resolution surface simulation at a certain hour takes 75 s. The dispersion kernel calculation using AERMET and AERMOD takes 3.4 s of this time. The preparation of the emission proxy data takes 6.6 s, mostly spent on interpolation and gridding of the traffic volumes. The surface concentration simulation takes 60.5 s, of which 96 % of the time is spent on Hadamard product calculations. Finally, the spatial interpolation of the surface measurements takes 4.2 s.

The emission optimisation is repeated once every simulation day. This takes 188 s if only surface measurements are considered (159 s are spent in emission proxy preparation). It takes ∼ 75 s longer if TROPOMI observations are also taken into account; time which is needed to perform the column simulations and the spatio-temporal interpolation to individual retrieval footprints.

In practice the computational time is less since the emission proxy calculations are shared between the simulation and the emission optimisation, requiring computation only once. The Hadamard product calculation scales linearly with the number of receptor points. For the surface concentration simulation this forms currently a computational bottleneck. However, as it involves straightforward matrix operations, the total computation time can be significantly reduced by parallelizing these tasks across multiple cores or GPUs.

4.3 Use of low-cost sensor data

As shown before in Mijling (2020), the Retina algorithm can also be applied to networks of low-cost sensors. Based on Bayesian principles in both emission optimisation and spatial assimilation, it effectively manages the greater inaccuracies associated with LCS data. The larger errors in the in-situ observations will slow down convergence to practical emission estimates in the optimisation phase (i.e., a longer lag period), but this is not necessarily a problem when emission trends are small over time.

Note that most NO₂ low-cost sensors used in current experimental networks suffer from creeping biases (Mijling et al., 2018; Li et al., 2021; WMO, 2023). Also reference measurements can be biased due to interfering gases (Steinbacher et al., 2007) or poor maintenance. Although always special care must be taken to remove this bias before application in Retina, integrating satellite measurements can help to reduce the introduced bias in the monitoring system.

5 Conclusion

The Retina algorithm has been designed to produce realistic high-spatiotemporal-resolution maps of urban air pollution based on heterogeneous air quality measurements. In this study, we implemented the updated Retina algorithm for NO₂ concentrations in Madrid and assessed the performance under different observation scenarios during March 2019. The updated algorithm, Retina v2, is faster and more accurate than its predecessor described in Mijling (2020). Most notably, it is now capable to incorporate tropospheric column concentrations retrieved with satellite instruments.

The use of proxy data for the description of urban emissions allows for convenient portability to other urban domains. Periodic emission optimisation guarantees that simulations match the observations, either satellite measurements, in-situ measurements, or both. Physics-based and running with modest computational power, Retina has comparable or better performance than data fusion methods depending on advanced, computational-demanding dispersion models, as well as machine learning approaches depending on extensive datasets.

When emissions are optimised using only TROPOMI measurements (representing the case of a city without an in-situ network), simulations of hourly NO₂ concentrations in March 2019 show a citywide RMSE of 19.3 µg m⁻³ with a bias of 0.8 µg m⁻³ . The algorithm's performance in March 2019 is broadly representative of annual averages, though seasonal variation exists. Correlation peaks at 0.86 in colder months with higher NO₂ levels, and drops to 0.70 in summer, when lower concentrations and reduced bias result in lower RMSE.

More accurate results are achieved when hourly in-situ measurements are used, as they allow for a better estimation of the diurnal emission cycle. However, if only a single station is available, and its measurements are biased or located in an area where dispersion modelling is problematic – due to e.g. incorrect proxy data – it can introduce systematic biases across the entire model domain. Incorporating satellite measurements or data from additional ground stations helps to reduce this bias.

The spatial interpolation of in-situ measurements in the simulation results improves the accuracy significantly: near observation sites it reduces simulation biases (e.g. due to inaccurate local emissions), and over larger distances it reduces simulation errors due to errors in background concentrations. Generally, including more stations leads to better results. Using all 24 ground stations in Madrid, the average correlation of hourly NO₂ time series increases to 0.90, with an RMSE of 13.0 µg m⁻³ , corresponding to a relative error of 36 %. Occasionally, the spatial interpolation introduces an extra bias. This suggests that there is further room for improvement in the covariance model used for interpolation.

The assimilation experiments for Madrid indicate that the added value of TROPOMI NO₂ measurements becomes negligible when hourly data from five or more ground-based stations at representative locations is available. However, this rule of thumb cannot be directly applied to other cities, as the contribution of TROPOMI depends on various factors, including city size, NO₂ concentration levels, and local meteorological conditions. Nevertheless, in many urban areas – especially those with sparse in situ monitoring – TROPOMI has the potential to provide substantial added value. Among approximately 2800 European cities with a population over 50 000, the European Environment Agency's AirBase (EEA, 2018) lists 2035 cities with at least one NO₂ monitoring station, but only 71 cities with five or more NO₂ stations (see Table S2).

The impact of satellite measurements in the Retina algorithm will be larger if observations at different times throughout the day could be included. Therefore, the next step will be preparation for data of the Sentinel-4 instrument aboard the MTG-S satellite, expected to be launched in late 2025. Operating from a geostationary orbit, Sentinel-4 will provide hourly measurements of NO₂ in Europe at 8 × 8 km² resolution with a revisit time of approximately 60 min. Once alternatives for CAMS background concentrations are available beyond the European domain, applications can extend to geostationary instruments such as GEMS (aboard GEO-KOMPSAT-2B, monitoring East Asia) and TEMPO (aboard Intelsat-40E, monitoring North America).