Improved constraints on ammonia emissions and deposition from co-assimilating NH₃ and NO₂ satellite observations over the Netherlands

Ammonia (NH₃) and nitrogen dioxide (NO₂) are key components of reactive nitrogen, strongly affecting air quality and ecosystem health. However, long-term constraints on ammonia emissions and deposition remain uncertain due to sparse in situ measurements and limitations of individual satellite products. We jointly assimilate five years (2018–2022) of NH₃ and NO₂ satellite observations over the Netherlands to improve constraints on reactive nitrogen concentrations, emissions, and deposition. NH₃ retrievals from the Infrared Atmospheric Sounding Interferometer (IASI) and the Cross-Track Infrared Sounder (CrIS) are combined with NO₂ observations from the TROPOspheric Monitoring Instrument (TROPOMI) within the LOTOS-EUROS chemical transport model using a Local Ensemble Transform Kalman Filter. The co-assimilation produces coherent year-to-year adjustments in modeled NH₃ concentration, emission, and deposition fields. Validation against measurements from the Dutch National Air Quality Monitoring Network (LML) shows reduced biases, clearer diurnal cycles, and improved correlations. Sensitivity experiments demonstrate that including TROPOMI NO₂ alongside IASI and CrIS NH₃ yields the lowest NH₃ surface bias vs. LML, highlighting the added value of coupling chemically related satellite observations. Comparisons with monthly Measurements of Ammonia in Nature (MAN) observations showed improved correlations but persistent spatial biases due to representativeness differences, while MAN sensors co-located with LML stations exhibited consistent improvements. These results demonstrate that co-assimilating complementary satellite observations can substantially improve constraints on ammonia emissions and deposition, with direct relevance for air-quality assessment and nitrogen policy applications.

1 Introduction

Ammonia (NH₃) and nitrogen dioxide (NO₂) are key reactive nitrogen species that play central roles in atmospheric chemistry and air quality. Through dry and wet deposition, these gases contribute substantially to the transfer of reactive nitrogen to the Earth's surface, where they influence both terrestrial and aquatic ecosystems. Although reactive nitrogen is essential for plant growth and ecosystem functioning, excessive deposition can cause environmental degradation, including soil acidification, eutrophication, and biodiversity loss (Galloway et al., 2008; Sutton et al., 2011).

NH₃ is the most abundant alkaline gas in the atmosphere, primarily emitted from agricultural activities such as livestock farming and fertilizer application (Erisman et al., 2011; Galloway et al., 2008). In agriculture, NH₃ is released through the decomposition of animal waste and the volatilization from nitrogen-rich fertilizers applied to soils (Erisman et al., 2008, 2011; Sutton et al., 2011; Paulot et al., 2014). This gas plays a vital role in the formation of secondary inorganic aerosols (SIAs) by reacting with acidic compounds such as sulfuric and nitric acids to form ammonium sulfate and ammonium nitrate (Erisman et al., 2011). SIAs contribute to fine particulate matter (PM_2.5), which affects atmospheric visibility and poses significant health risks, including respiratory and cardiovascular diseases (Pope and Dockery, 2006; Pope et al., 2009).

NO₂ in the atmosphere originates primarily from nitric oxide (NO), which is emitted during fossil fuel combustion in vehicles, power plants, and industrial processes and rapidly oxidized in air. Together, NO and NO₂ form the NO_x family, a key precursor to tropospheric ozone and nitrate aerosols (Crutzen, 1979; Ojeda-Castillo et al., 2025; Zara et al., 2021). NO₂ participates in photochemical reactions leading to the formation of ozone and particulate nitrates, contributing to air pollution and smog formation (Zhang et al., 2011). In addition, NO₂ indirectly influences the atmospheric lifetime of NH₃ by modulating the production of nitric acid (HNO₃): increased oxidation of NO₂ enhances HNO₃ formation, which in turn promotes the partitioning of NH₃ into ammonium nitrate. Elevated levels of NO₂ and its reaction products are associated with respiratory problems such as asthma and decreased lung function, and contribute to environmental degradation through acid deposition and nutrient imbalances in ecosystems (Eum et al., 2022; Jonson et al., 2017; Liu et al., 2021).

Understanding the spatial and temporal distributions of NH₃ and NO₂ is essential for accurate atmospheric modeling and the development of effective emissions control strategies and regulations. However, the reactive nature and short atmospheric lifetimes of NH₃ and NO₂ pose significant challenges to monitoring and modeling of these species. Emissions inventories derived from bottom-up approaches often have large uncertainties associated with them, particularly in the case of NH₃ where these uncertainties can be greater than 100 % for some emissions sectors (Kuenen et al., 2022). Ground-based measurement networks provide high-accuracy observations but differ greatly in spatial density; NO₂ is routinely monitored at hundreds of sites across Europe, whereas NH₃ is measured at far fewer locations, resulting in much sparser spatial coverage and limited temporal resolution. Moreover, these networks generally measure surface concentrations and do not provide full atmospheric column information. Space-based observations from instruments such as IASI on the MetOp series, CrIS on the Suomi-NPP and NOAA-20/21 platforms, and TROPOMI on Sentinel-5 Precursor offer broader spatial coverage, but with limited temporal sampling due to discrete overpass times. Satellite measurements of this type can be subject to errors resulting from, for example, instrumental noise, retrieval biases, and residual cloud cover.

Inverse modeling methods present a powerful approach to overcome these limitations by integrating observational data into chemical transport models to yield optimized estimates of atmospheric constituents. By assimilating large numbers of NH₃ and NO₂ measurements, it is possible to improve the accuracy of simulated concentrations, diurnal cycles, and associated emission and deposition fields. Assimilating NO₂ observations in particular helps to better constrain the NO_x budget, which in turn reduces uncertainties in the formation and loss of nitric acid (HNO₃) and thus indirectly constrains the NH₃ fields through their chemical coupling. Through such an approach, air quality predictions can be improved and policy-relevant insights can be gained for mitigating the impacts of reactive nitrogen on ecosystems and human health.

In this paper, we perform a co-assimilation of measurements of NH₃ from IASI, CrIS and NO₂ from TROPOMI in the LOTOS-EUROS local ensemble transform Kalman filter (LETKF) over a model domain encompassing the Netherlands and adjacent parts of northwestern Germany. This is a particularly relevant region for studying atmospheric NH₃, as it forms a major reactive nitrogen hot-spot in Europe due to intensive agriculture, especially livestock production and fertilizer use. At the same time, accurate quantification of NH₃ is particularly important for the Netherlands given the ongoing nitrogen crisis and the associated pressures of nitrogen deposition on sensitive ecosystems. Including the neighboring German source regions is also necessary because NH₃ and secondary inorganic nitrogen are influenced by cross-border transport, such that concentrations and deposition over the Netherlands cannot be interpreted from domestic emissions alone. We evaluate resulting optimized emissions and deposition fields, with a focus on NH₃, and we compare the results against independent observations from ground-based measurement networks.

2 Methodology

2.1 The LOTOS-EUROS model

In this study, we utilized the LOTOS-EUROS (LOng Term Ozone Simulation-EURopean Operational Smog) v2.2.009 chemical transport model to simulate atmospheric concentrations over the study region (Manders et al., 2017). LOTOS-EUROS is a three-dimensional Eulerian model designed for regional air quality assessments and operational forecasting in Europe. It effectively simulates the dispersion, chemical transformation, and deposition of atmospheric pollutants, including gases and aerosols. LOTOS-EUROS is part of the Copernicus Atmospheric Monitoring Service (CAMS) European air quality ensemble (Colette et al., 2025). This service provides forecasts for the main air pollutants using an ensemble of state-of-the-science CTMs. Within CAMS, LOTOS-EUROS is regularly validated against in-situ observations and TROPOMI satellite data, as well as evaluation against the other ensemble members (Peuch et al., 2022; Colette et al., 2025). LOTOS-EUROS also has participated in numerous model inter-comparisons, typically showing a strong performance (Bessagnet et al., 2016; Colette et al., 2017; Vivanco et al., 2017).

The model incorporates detailed representations of atmospheric processes such as advection, diffusion, bi-directional fluxes, and chemical reactions. It uses the Carbon Bond Mechanism IV (CBM-IV) for gas-phase chemistry, which includes a comprehensive set of reactions relevant to ozone formation and other photochemical oxidants. Aerosol dynamics are modeled using size-resolved modules that account for primary emissions, secondary formation, and processes like coagulation and deposition.

Model outputs, including concentrations of key pollutants such as NH₃, NO₂, ozone (O₃), and particulate matter (PM₁₀ and PM_2.5), are regularly validated against observational data from the Dutch LML air quality monitoring network, measurements from the German environmental agency (the Umweltbundesamt; UBA) and the EBAS network throughout the European Union. The NH₃, NO₂ and particulate matter components from the model are also frequently validated against the EMEP model and other models within the CAMS model ensemble (Schaap et al., 2008; Manders et al., 2017; Tsyro et al., 2022).

2.1.1 The Local ensemble transform Kalman filter

The Ensemble Kalman filter (EnKF; Evensen, 2003) is a sequential data assimilation method in which uncertainties in the model state are represented by an ensemble of simulations, which is updated using observations. In this study, we use the Local Ensemble Transform Kalman Filter (LETKF), a localized formulation of the EnKF that updates the model state by combining the ensemble forecast with observations within a defined spatial neighborhood. The LOTOS-EUROS LETKF v3.0.7 applied here has previously been used in studies of particulate matter (Lopez-Restrepo et al., 2020), NO₂ (Timmermans et al., 2019), and NH₃ (Van Der Graaf et al., 2022). The formulation is described in Hunt et al. (2007), and the implementation used here follows Shin et al. (2016). Compared with the standard Ensemble Kalman Filter, the LETKF is more computationally efficient and performs the analysis on a per-grid-cell basis using only nearby observations, which are determined using a specified spatial localization length (discussed in further detail in Sect. 2.1.2).

In the present application, the LETKF uses an augmented state vector,

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\mathbf{X}=(\mathit{c},\mathit{\beta }),\end{array}$$

where c denotes the three-dimensional trace-gas concentration fields and β denotes one or more two-dimensional emission perturbation fields for the optimized species. The primary objective of the filter is to estimate β, which defines multiplicative per-grid-cell scaling factors applied to the prior emissions. The perturbed emissions are computed as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\mathit{E}=max\left(\mathrm{0},{\mathit{E}}_{\text{base}}\cdot \mathit{\beta }\right),\end{array}$$

where E_base is the prior emission field. This clipping ensures that the updated emissions remain non-negative and prevents non-physical emission sinks. The concentration field c is included in the augmented state because it provides the dynamical link between the emission perturbations and the observations. During the forecast step, the model propagates the concentration fields forward using emissions modified by the current β values, so that free-running simulations, or simulations in regions without recently assimilated observations, remain influenced by earlier analysis updates.

The temporal variability in the emissions is specified in the assimilation through the temporal correlation length τ (set to 3 d for NH₃ based on Van Der Graaf et al., 2022, and 1 d for NO₂), with the temporal correlation coefficient α defined as Lopez-Restrepo et al. (2020), Van Der Graaf et al. (2022):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}_k={\mathrm{e}}^{-|t_k-t_{k-\mathrm{1}}|/\mathit{\tau }}\end{array}$$

where t_k and t_k−1 are successive hourly analysis times. This formulation ensures that, in the absence of new observations for an extended period, the influence of past updates diminishes and the system progressively returns toward the a priori emission state. The LETKF analysis is applied at hourly analysis times throughout the simulation period. At each analysis time t_k, only the emission state corresponding to the current time step is updated; emissions from previous time steps are not retrospectively adjusted. The temporal correlation coefficient α_k in Eq. (3) therefore does not define a multi-time assimilation window, but instead controls the persistence of emission adjustments between successive hourly analysis times.

In the main co-assimilation configuration used in this study, species-specific emission perturbation factors are optimized simultaneously for both NH₃ and NO₂. As a result, assimilated NO₂ observations constrain the NO_x emission field directly rather than allowing their signal to be attributed to NH₃ emissions alone. The uncertainty in this system is represented by an ensemble of N members. In this study, N=12 was used following Van Der Graaf et al. (2022), who applied the same LOTOS-EUROS LETKF framework. This relatively modest ensemble size is feasible here because the analysis is strongly localized in space and the short atmospheric lifetimes of NH₃ and NO₂ lead to comparatively compact, local covariance structures. For the initialization of the ensemble, the emission perturbation factors are sampled from a normal distribution with a mean of 1.0 and a standard deviation of 0.5, corresponding to 50 % uncertainty around the prior emissions. The resulting emissions are then constrained through Eq. (2), such that perturbation factors that would otherwise produce negative emissions instead yield zero.

The LETKF operates in two sequential steps, the forecast and the analysis. In the forecast step, the state vector ensemble ${\mathbf{X}}^f=\mathit{\{}{\mathit{x}}_{\mathrm{1}}^f,{\mathit{x}}_{\mathrm{2}}^f,\mathrm{\dots },{\mathit{x}}_N^f\mathit{\}}$ is propagated forward in time using the model dynamics, where ${\mathit{x}}_i^f$ represents the ith forecast ensemble member. Each analyzed ensemble member at time t_k−1 is propagated forward to t_k according to Van Der Graaf et al. (2022):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{\mathit{x}}_i^f\left(t_k\right)=M_{k-\mathrm{1}}\left({\mathit{x}}_i^{\mathrm{a}}\left(t_{k-\mathrm{1}}\right)\right),\end{array}$$

where the model operator M_k−1 describes the forward model simulation from t_k−1 to t_k, including the application of the emission perturbation factors to the prior emissions and the persistence of β between successive analysis times through the temporal correlation coefficient α_k. Here, ${\mathit{x}}_i^{\mathrm{a}}\left(t_{k-\mathrm{1}}\right)$ denotes the ith analyzed ensemble member at time t_k−1. Once the ensemble has been propagated forward in time, the ensemble mean ${\stackrel{\mathrm{‾}}{\mathit{x}}}^f$ and forecast error covariance P^f are calculated as

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{x}}}^f={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\mathit{x}}_i^f,\text{(6)}& {\displaystyle }& {\displaystyle }{\mathbf{P}}^f={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{i=\mathrm{1}}^N\left({\mathit{x}}_i^f-{\stackrel{\mathrm{‾}}{\mathit{x}}}^f\right){\left({\mathit{x}}_i^f-{\stackrel{\mathrm{‾}}{\mathit{x}}}^f\right)}^{\mathrm{T}},\end{array}$$

When new observations y^obs are available, the analysis ensemble X^a is obtained by updating each ensemble member according to Van Der Graaf et al. (2022):

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{\mathit{x}}_i^{\mathrm{a}}={\mathit{x}}_i^f+{\mathbf{P}}^{\mathrm{a}}{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\left({\mathit{y}}^{\text{obs}}-\mathit{h}\left({\mathit{x}}_i^f\right)\right).\end{array}$$

Here, $\mathit{h}\left({\mathit{x}}_i^f\right)$ denotes the model-simulated equivalent of the satellite retrieval, H is the linearized observation operator, R is the observation error covariance matrix, and P^a is the analysis error covariance, which is computed from Shin et al. (2016):

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathbf{P}}^{\mathrm{a}}={\left[{\mathbf{P}}^f{\mathbf{H}}^{\mathrm{T}}{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}+\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{P}}^f.\end{array}$$

The simulated satellite observations are computed following the averaging kernel formalism of Rodgers (2000) and Rodgers and Connor (2003). First, the model state is interpolated to the retrieval grid using the gridding operator G:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\mathit{x}}_{r,i}=\mathbf{G}{\mathit{x}}_i^f.\end{array}$$

To ensure that the comparison with the satellite product is made at the same effective vertical resolution as the retrieval, the averaging kernel A is applied:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }\mathit{h}\left({\mathit{x}}_i^f\right)={\mathit{x}}_a+\mathbf{A}\left({\mathit{x}}_{r,i}-{\mathit{x}}_a\right),\end{array}$$

where x_a is the a priori profile used in the retrieval. In the linear approximation, the observation operator becomes:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\mathbf{H}=\mathbf{AG}.\end{array}$$

The observations are assumed to satisfy

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathit{y}}^{\text{obs}}=\mathit{h}\left({\mathit{x}}^{\text{true}}\right)+\mathit{v},\mathit{v}\sim \mathcal{N}(\mathrm{0},\mathbf{R}),\end{array}$$

with R taken from the retrieval error covariance matrices of the satellite data products and represents both measurement and representativeness errors (Van Der Graaf et al., 2022). Once the update is complete, the analyzed ensemble X^a becomes the initial condition for the next forecast step.

2.1.2 Spatial localization

To ensure computational efficiency and avoid spurious correlations, the LETKF applies spatial and temporal localization in a per-grid-cell approach following Shin et al. (2016). In contrast to approaches that apply covariance localization directly to the background error covariance matrix, the LETKF implementation used here applies localization in observation space by selecting and weighting nearby observations for each local analysis. As a result, the localization length in this framework should be interpreted as application-specific and is not expected to match the much larger values used in numerical weather prediction studies such as Shin et al. (2016), where the analyzed variables exhibit broader synoptic-scale spatial correlations. The temporal localization was described in the previous section, and is applied using Eq. (3). For the spatial localization, the simulated observations are first computed for all ensemble members, and then for the given grid cell to be analyzed, all observations (both simulated and real) within 3.5ρ distance are gathered, with ρ being the specified localization length (in units of km). The analysis is then performed using the collected observations, with the contribution of each observation to the local analysis decreasing smoothly with distance from the analyzed grid cell. In the present implementation, this distance weighting is represented using a Gaussian decay function:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }w\left(\mathrm{\Delta }d\right)=\exp \left(-{\displaystyle \frac{\mathrm{\Delta }{d}^{\mathrm{2}}}{\mathrm{2}{\mathit{\rho }}^{\mathrm{2}}}}\right),\end{array}$$

where Δd is the distance (in km) between the observation and the model grid point. Thus, observations closest to the analyzed grid cell have the largest influence, while the contribution of more distant observations decreases smoothly with distance. This is consistent with the localized observation-space LETKF framework described by Shin et al. (2016), although the functional form shown here corresponds to the implementation used in the present study. In this study, the spatial correlation lengths are chosen to be consistent with the horizontal representativeness of the corresponding satellite retrievals and are therefore guided by the mean footprint sizes. This also reflects the relatively local emission-concentration relationships of short-lived reactive gases such as NH₃ and NO₂. Using a localization length much larger than the footprint would allow a single observation to influence the analysis over spatial scales that are not resolved by the measurement, potentially producing spurious long-range increments and unrealistically smooth updates. To ensure spatial consistency between the retrieval resolution and its influence in the LETKF, we use ρ=15 km for CrIS and IASI NH₃ and a smaller ρ=5 km for TROPOMI NO₂, reflecting the higher spatial resolution of the latter.

2.1.3 Model configuration

In this section, a short summary of the most important model inputs and configuration parameters are provided. For a more detailed description of the LOTOS-EUROS model we direct the reader to Manders et al. (2017). The LETKF v3.0.7 is coupled to version 2.2.009 of the LOTOS-EUROS model. An initial long simulation covering the period of 2018–2022 was performed on a domain that covers the majority of Europe (15° W–35° E; 35–70° N) with a resolution of approximately 25 km×25 km, and the output from this run was utilized as boundary conditions for the assimilation run which was performed on a domain covering most of North-western Europe (2–16° E; 47–56° N) with a resolution of 7 km×7 km. The model is driven by meteorological fields obtained from the ECMWF short-term forecast model at a 3-hourly temporal resolution which is then interpolated to an hourly frequency within the model. The simulations were conducted using 12 vertical levels, extending from the ground to about 10 km above the earth's surface, matching the vertical layer structure of the ECMWF meteorology dataset.

The base emission dataset used in the model simulations combines the European-scale CAMS-REG-v5.1 inventory (Kuenen et al., 2022) with higher-resolution national inventories for the Netherlands (Emissieregistratie; ER) and Germany (Gridding Emission Tool for ArcGIS; GRETA). This combined CAMS GRETA-ER emissions dataset was developed within the National Kennisprogramma Stikstof (NKS) funded by the Dutch Ministry of Agriculture, Fisheries, Food Security and Nature (LVVN). The base emission dataset was compiled using the corresponding inventory year where available; for years after 2019, the 2019 emission totals were used as the baseline because 2019 was the most recent year available in the harmonized CAMS-GrETa-ER emissions dataset, and these emissions were then adjusted dynamically according to meteorological conditions. As a result, year-to-year variations in the base emissions remain relatively small and primarily reflect meteorological influences rather than structural changes in activity or policy. Temporally, the emissions are distributed using hourly time factors specific to aggregated source categories. For agricultural NH₃ emissions, a meteorologically dependent parameterization is applied that accounts for weather-driven shifts in fertilizer application timing, following the approach of Ge et al. (2020). Emissions are also vertically distributed according to sector-specific release heights, which is particularly relevant for industrial and power-generation sources where average stack heights determine the initial plume elevation.

The model output for all major nitrogen species consists of simulated concentrations (surface, and all vertical layers) and wet and dry deposition, as well as concentrations matched at the footprints of all satellite products (i.e. IASI, CrIS, TROPOMI). To match the model with the satellite footprints, and for the ingestion of these observations in the LETKF, the CAMS Satellite Operator (CSO) is used. CSO (https://ci.tno.nl/gitlab/cams/cso, last access: 6 July 2026) is an open-access tool developed at the Netherlands Organisation for Applied Scientific Research (TNO) and implemented to facilitate fast intercomparisons between modelled and satellite concentrations. The tool consists of two entities: a pre-processor to download, select, and convert satellite observations into a common format, accompanied by post-processing tools to aggregate and visualize the data; and a source code that can be used within regional air quality modelling and assimilation systems such as LOTOS-EUROS. The CSO module is able to read the files created by the pre-processor, simulate satellite observations using model variables, and apply observational operators where applicable.

2.2 Satellite datasets

In the LETKF configuration used in this study, the assimilated satellite observations consist of NH₃ total columns from IASI and CrIS, and NO₂ tropospheric vertical column densities from TROPOMI.

2.2.1 IASI NH₃

The IASI instruments onboard the MetOp-A, -B, and -C satellites are in sun-synchronous orbits, with Equator crossing times at approximately 09:30 and 21:30 LT (Clerbaux et al., 2009). Although the three platforms fly in the same local-time orbit, they are phased along that orbit and therefore do not acquire measurements simultaneously over the same ground location; the temporal separation between the platforms is on the order of 45 min. The data products of IASI-A, -B, and -C span the periods October 2007–October 2021 (IASI-A), March 2013 onward (IASI-B), and September 2019 onward (IASI-C), with the latter two instruments still operational. Each of the three IASI instruments has an observational swath width exceeding 2000 km, with a pixel footprint of approximately 12 km in diameter at nadir, increasing to around $\sim \mathrm{20}\mathrm{km}\times \mathrm{40}\mathrm{km}$ at the swath edges. In this study, we utilize the most recent IASI product, the Artificial Neural Network for IASI (ANNI)v4 (Clarisse et al., 2023), which is an updated version of the earlier IASI-NNv2 and LUT products (Van Damme et al., 2014, 2017).

Similar to its predecessors, the IASI-ANNIv4 retrieval involves two steps. First, the hyperspectral range index (HRI) is calculated to characterize the NH₃ signal strength in each spectrum (Van Damme et al., 2017; Clarisse et al., 2023). The second step utilizes a neural network trained on a large dataset of modeled data, which links the HRI to the NH₃ total column. The primary improvement in this version is the addition of a column averaging kernel to facilitate comparisons with in-situ and model data, enabling the effects of the a-priori profile shape to be considered. Earlier ANNI NH₃ products were evaluated in previous studies by Dammers et al. (2016), Guo et al. (2021), Kutzner et al. (2021), Herrera et al. (2022), including against ground-based FTIR observations from the Network for the Detection of Atmospheric Composition Change (NDACC). The validation of IASI NH₃ against NDACC measurements by Dammers et al. (2016) showed on average favorable correlations (R=0.80), with a mean low bias on the order of 35 %.

In the current study, we apply the recommended data quality filters, namely, with the pre-filter and post-filter set to 1, and only observations with a cloud_fraction <25 % are used.

2.2.2 CrIS NH₃

The CrIS-1 instrument on the Suomi-NPP and CrIS-2 on NOAA-20 were launched in October 2011 and November 2017, respectively, with CrIS-Fast Physical Retrieval (CFPR) NH₃ data products starting in May 2012 and March 2019. Both instruments are in sun-synchronous orbits, providing global coverage twice daily at around 13:30 and 01:30 local solar time, with overpasses within 45 min of each other. They offer observations with circular pixel footprints of approximately 14 km at nadir over a 2200 km wide swath. This study utilizes the CFPR NH₃ product version 1.6.4 (Shephard and Cady-Pereira, 2015; Shephard et al., 2020). The CFPR method involves a physical retrieval based on the Rodgers (2000) optimal estimation method combined with a fast optimal spectral sampling forward model (Moncet et al., 2008), minimizing the residual between measured and simulated spectra. The excellent signal-to-noise ratio (∼1600) of the CrIS instrument in the ammonia spectral region enables detection sensitivities of approximately 0.5 ppbv near the surface, or 3.5×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$ for total columns, under typical atmospheric conditions (∼50 % detection rate). Under highly favorable infrared remote sensing conditions (e.g. strong thermal contrast) the detection limit can improve to about 0.2 ppbv, corresponding to a ∼10 % detection rate (Shephard et al., 2025). Since the last major validation study by Dammers et al. (2017), the product has undergone several iterations, including the addition of a cloud flag based on VIIRS data, non-detects, and a quality flag (White et al., 2023). The validation study by Dammers et al. (2017) reported a good correlation between FTIR and satellite observations (R≈0.8) with a slight high bias (slope=1.02). For higher column concentrations, CrIS observations showed a small positive difference with the ground-based FTIR measurements around 25 %–50 %, while for lower concentrations, the bias increased to 2.5×10¹⁵ $\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$ with a standard deviation of around 50 %–100 %. This study only includes observations with a quality_flag of ≥3, thereby excluding failed or lower-confidence retrievals, and with a cloud_flag equal to 0 (clear-sky scenes). Only daytime observations are assimilated. It should be noted that the CrIS-1 instrument suffered a failure of its mid-wave IR (MWIR) band from 26 March to 24 June 2019, leading to a data gap in this period (Iturbide-Sanchez et al., 2022).

2.2.3 TROPOMI NO₂

The TROPOMI instrument, aboard the Sentinel-5 Precursor (S5P) polar-orbiting satellite, is a nadir-viewing spectrometer designed for atmospheric observations. It crosses the equator at approximately 13:30 LT. The instrument measures radiation across the ultraviolet, visible, and infrared spectral ranges, enabling the monitoring of atmospheric trace gases and aerosols (Veefkind et al., 2012). TROPOMI has a swath width of approximately 2600 km, and the NO₂ product has a nadir spatial resolution of 7.2 km in the along-track direction and 3.6 km in the across-track direction, improving to 5.6 km×3.6 km after 6 August 2019 (van Geffen et al., 2022). The retrieval of NO₂ columns follows a three-step process. First, the NO₂ slant column density is calculated from the L1b spectra recorded by TROPOMI using a DOAS fitting algorithm. This slant column is then separated into stratospheric and tropospheric components through data assimilation using the TM5-MP model, which operates at a horizontal resolution of 1°×1° (Williams et al., 2017). Finally, slant column densities are converted to vertical column densities (VCD) by applying total and altitude-dependent air mass factors (AMFs). These AMFs are influenced by several factors, including NO₂ vertical profiles obtained from TM5-MP, the satellite's viewing geometry, surface albedo, surface pressure, and characteristics of clouds and aerosols. Further details on the retrieval process can be found in van Geffen et al. (2022) and in the algorithm theoretical baseline document (van Geffen et al., 2024).

Routine validation of TROPOMI NO₂ measurements against ground-based MAX-DOAS observations from 29 stations has revealed a mean bias of −28 %, increasing to −40 % in regions with heavy pollution (Lambert et al., 2024). This bias largely stems from the TM5-MP vertical profiles, which inadequately resolve high-concentration hot-spots and show deviations in the profile shape, particularly near the surface (Chan et al., 2020; Verhoelst et al., 2021). To address these discrepancies, the a-priori vertical profile can be updated using one derived from a higher-resolution air quality model, which has been shown to partially mitigate biases (Griffin et al., 2019; Zhao et al., 2020; Judd et al., 2020; Douros et al., 2023). This correction process employs TROPOMI averaging kernels and is described in detail in the TROPOMI NO₂ Product User Manual (Eskes et al., 2024).

In this study, we used the VCDs from the reprocessed TROPOMI NO₂ version 2.4.0 dataset. To ensure data reliability, observations with a quality assurance value below 0.75 were excluded. This threshold effectively eliminates pixels with cloud radiance fractions exceeding 0.5, thereby reducing the impact of uncertain retrievals (van Geffen et al., 2022).

2.3 Ground-based in-situ measurements

2.3.1 The LML network

The Dutch National Air Quality Monitoring Network, known as the “Landelijk Meetnet Luchtkwaliteit” (LML) (Elskamp, 1989; Elzakker and Buijsman, 1999), is a comprehensive ground-based measurement network designed to monitor air quality across the Netherlands. Operated by Rijksinstituut voor Volksgezondheid en Milieu (RIVM), the LML network consists of a large number of monitoring stations distributed across urban, suburban, and rural areas. These stations continuously collect data on various air pollutants, including PM₁₀, PM_2.5, NO₂, NH₃, O₃, sulfur dioxide (SO₂), carbon monoxide (CO), and volatile organic compounds (VOCs).

The network provides real-time data, which is crucial for assessing the air quality in different regions and understanding the impact of pollution on public health and the environment. LML stations employ state-of-the-art sensors and analytical techniques to ensure high data accuracy and consistency, enabling authorities to monitor trends, detect exceedances of air quality standards, and develop policy interventions when necessary. NH₃ measurements are made using the miniDOAS, an active instrument that utilizes the differential optical absorption spectroscopy (DOAS) measurement technique. These instruments have an NH₃ detection limit of roughly 0.25 µg m⁻³ and an estimated precision of 0.1 µg m⁻³ for hourly averaged observations (Berkhout et al., 2017), however, no measurement uncertainties are provided in the dataset.

Data from the LML are publicly accessible (http://www.luchtmeetnet.nl, last access: 6 July 2026), allowing citizens, researchers, and policymakers to track air quality levels in near real-time. This transparency helps raise awareness about air pollution issues and supports efforts toward improving air quality across the country. The LML network is also integrated with broader European air quality initiatives, contributing to the wider understanding of transboundary pollution and climate change mitigation efforts. The LML NH₃ measurements were previously applied to study long-term trends in the Netherlands by van Zanten et al. (2017).

Ten LML sites were selected in total for the comparisons with the model simulations, and the coordinates, details, and species measured for each of these sites are provided in Table 1. The locations of the sites within the Netherlands are shown on a map in Fig. 1. Six of the chosen LML sites provide hourly NH₃ surface concentration measurements; De Zilk-Vogelaarsdreef, Valthermond-Noorderdiep, Vredepeel-Vredeweg, Wekerom-Riemterdijk, Wieringerwerf-Medeblikkerweg, and Zegveld-Oude Meije. These sites were selected because they provided hourly NH₃ measurements over a long period (>1 year). In addition to the NH₃ measurements, wet deposition measurements of dissolved ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) concentrations in precipitation are made periodically (i.e. at irregular intervals) at Biest Houtakker-Biestsestraat, De Bilt-Wilheminalaan, Philippine-Stelleweg, Speuld-Garderenseweg, Valthermond-Noorderdiep, Vredepeel-Vredeweg, and Wieringerwerf-Medemblikkerweg. Wet deposition is also measured at the De Zilk-Vogelaarsdreef site, which serves as a European Monitoring and Evaluation Programme (EMEP) location where dissolved ammonium is monitored on a daily basis. The observed ${\mathrm{NH}}_{\mathrm{4}}^{+}$ measurements were paired with the model output and converted to monthly mean fluxes using the corresponding measured and modeled precipitation amounts. The modeled and observed precipitation agreed well on average, though some transient mismatches occurred. To minimize the influence of transient precipitation mismatches on the wet-deposition comparison, a relatively strict filter was applied: cases where the mean absolute deviation of the measured and modeled precipitation differed by more than 1σ were excluded, as were measurements with very low precipitation (<0.1 mm). Sensitivity tests with looser thresholds (2σ and 3σ) led to poorer agreement and increased spread in the deposition comparison, so the 1σ filter was adopted for the final wet-deposition evaluation.

Table 1 Locations and details of LML sites used for comparisons with the LOTOS-EUROS LETKF simulations.

Figure 1 Map of the locations of the selected LML sites in the Netherlands used for comparisons with the model simulations. Underlying basemap data sourced from OpenStreetMap contributors (2017).

The NH₃ measurement sites occasionally experience local-scale enhancement events (i.e. nearby manure dumping events), which are not spatially representative of the model resolution. For 2020 onward, official flags are provided in the LML dataset to enable the filtering of such events, but for the years prior to this no such information is available. As a result, additional filtering was manually applied, and any hourly NH₃ observations with concentrations exceeding 100 µg m⁻³ are removed and excluded from the comparisons.

2.3.2 The MAN network

The “Measuring Ammonia in Nature” (MAN) network was established in 2005 to monitor atmospheric ammonia concentrations in nature reserve areas across the Netherlands, with a particular focus on nitrogen-sensitive Natura 2000 areas (Lolkema et al., 2015). The network provides essential data for assessing national ammonia concentration trends, validating air quality models, and analyzing regional variability.

The network employs commercial Gradko passive ammonia samplers, which are cost-effective, easy to deploy, and well-suited for large-scale monitoring (Lolkema et al., 2015; Noordijk et al., 2020). These samplers are calibrated monthly using active ammonia sampling devices within the LML network (Lolkema et al., 2015). Local volunteers, often conservation wardens, handle the monthly exchange of samplers, ensuring consistent data collection across diverse habitats.

Currently, the MAN network includes over 300 sampling sites across the Netherlands. The ammonia concentration data gathered by the network facilitates the identification of spatial concentration patterns and regional anomalies throughout the country. Note that in comparison to LML, the MAN samplers are typically placed in nature areas and away from source regions.

The MAN network's data plays a critical role in assessing the effectiveness of environmental policies aimed at reducing nitrogen emissions and in analyzing long-term trends in nitrogen deposition. The network can detect annual trends as low as 3 % over extended time series, making it a valuable tool for air quality management and biodiversity conservation in nitrogen-sensitive areas (Lolkema et al., 2015).

To reduce meteorological influences on the passive samplers, the MAN network is calibrated against high-performance reference instruments for NH₃ at six locations in the Netherlands (see Berkhout et al., 2017). The calibration process is detailed in Noordijk et al. (2020). The uncertainty in the MAN measurements consists of two components: the measurement uncertainty and the calibration uncertainty. The random uncertainty is approximately 0.9 µg m⁻³ for a single monthly value and 0.32 µg m⁻³ for annual averages. The systematic uncertainty is estimated at 28 % for monthly values and 10 % for annual averages (Noordijk et al., 2020).

For comparison with the LETKF simulations, all available MAN data from the period 2018–2022 was used, comprising a total of 309 standard sites, and 6 additional MAN calibration sensors located at LML measurement sites. While not all sites provide uninterrupted time series over the full period, the large number of sites and monthly measurements ensures statistically robust results.

3 Results and Discussion

3.1 Optimized emission fields

We first examine the NH₃ emission fields pre- and post-assimilation to evaluate the impact of ingesting satellite observations on the model simulation. The base and LETKF-optimized NH₃ yearly total emissions for each individual year and for the mean over 2018–2022 are shown in Fig. 2. Unless otherwise stated, the LETKF-optimized simulation refers to the main co-assimilation run using NH₃ observations from IASI and CrIS together with NO₂ observations from TROPOMI. The relative-difference plots in Fig. 2c reveal a consistent spatial pattern: increases across much of the south and east of the Netherlands and decreases in the north and the adjacent regions of Germany. The largest mean increase occurred in 2020 (+14.3 %), the smallest in 2021 (+3.0 %), with a period-average change of +8.0 %. The relative-difference panels in Fig. 2c should therefore be interpreted together with the underlying base and optimized emission fields in Fig. 2a and b, since large percentage changes can still correspond to modest absolute changes where baseline emissions are low. It should be noted that 2020 was an exceptional year, with unusually warm and sunny summer months, which in turn led to increased volatilization of NH₃ and higher emissions. This is likely one of the key factors contributing to the higher emissions change post-assimilation in that particular year.

Figure 2 (a) The total NH₃ emissions by year and the mean emissions 2018–2022 period from the base CAMS GRETA-ER inventory, (b) the same but for the optimized emissions from the LETKF analysis, and (c) the mean relative differences between the base and the optimized emission fields.

The persistent increases in the south-eastern Netherlands mirror the patterns reported by Ge et al. (2020), who found that incorporating detailed agricultural activity data into the Monitoring Atmospheric Composition and Climate (MACC) inventory led to higher emissions in this region due to improved spatial allocation of sources, more accurate representation of manure management and application timing, and region-specific regulatory constraints. They also showed that emissions here are disproportionately influenced by intensive pig and poultry farming, whereas dairy cattle dominate in much of the rest of the country. In contrast, the decreases we find in the north are consistent with their observation that refined allocation can reduce emissions where earlier inventories overestimated activity, particularly for dairy cattle. These parallels suggest that the positive adjustments in our assimilation likely reflect structural biases in the base inventory, both in the spatial distribution and the livestock-sector partitioning of emissions, rather than being solely an artifact of the assimilation.

Meteorological effects can also have a significant impact on NH₃ emissions, with temperature, precipitation, and wind speed strongly influencing volatilization rates and atmospheric transport (Søgaard et al., 2002; Gyldenkærne et al., 2005; Ge et al., 2023). Warm, dry conditions can substantially enhance emissions above climatological norms, while precipitation events may suppress them or trigger post-event peaks. Such variability may contribute to the year-to-year differences in our assimilation adjustments, and could help explain the particularly large increase in 2020, when the Netherlands experienced a very sunny spring and an unusually warm, dry summer.

Figure 3 A time-series of monthly total NH₃ emissions calculated over the model domain shown in Fig. 2 for the 2018–2022 period from the base simulation and the optimized LETKF simulation. The relative difference (optimized – base)/base in % is shown by the dashed line.

In addition to the annual emission changes, it is also informative to examine the temporal evolution of emissions at finer timescales. Figure 3 presents the time series of monthly NH₃ emission totals, aggregated over the same region shown in Fig. 2, and provides further insight into how the assimilation influences variability across individual months. The base emissions show a similar seasonal cycle in all years, characterized by a pronounced spring peak and a smaller secondary peak in summer. This double-peaked seasonal cycle likely reflects the combination of the prescribed seasonal timing in the agricultural NH₃ emission parameterization and meteorologically driven variability in NH₃ volatilization under warmer conditions. Relative to the base simulation, the optimized emissions show a broadly consistent pattern of changes across all years, with a reduction in the spring peak and a corresponding increase during the summer months. Particularly large decreases in the spring peaks are seen for 2018–2020, whereas only a minimal decrease occurs in 2021 and a slight increase is found in 2022. As seen in Fig. 3, 2020 exhibits the largest emissions changes, with particularly large increases in the emissions (on the order of +70 %) in the LETKF-optimized simulation between April and September. In contrast, 2021 shows the smallest emissions changes, with a smaller-than-average increase in the emissions in the optimized simulation in May and June, but mostly minor changes throughout the remainder of the year.

Similar NH₃ seasonal cycles to that of the optimized simulation have been presented in recent modeling and observational studies of the region. The updated temporal emission profiles presented by Ge et al. (2020), and which are implemented in the model simulations used in this study, show shifts in the early-year emission peaks to later in the springtime, and increases in the overall summertime emissions. A long-term time-series (2008–2020) of IASI NH₃ total column measurements over the Netherlands in Van Damme et al. (2022) displays a springtime peak in April and a slightly smaller but somewhat comparable secondary peak in July to August. A similar pattern of reduced spring peaks and enhanced summertime emissions was also found by Ding et al. (2024) after assimilating CrIS NH₃ observations into the DECSO system. In addition, a recent two-year field campaign in key agricultural regions of the Netherlands by Lô et al. (2025) observed strong peaks in surface NH₃ concentrations during the summer months, driven by increased volatilization under warmer conditions. Despite the inclusion of the Ge et al. (2020) temporal profiles, the model simulations in this study still underestimate summertime emissions relative to these independent datasets, suggesting that the seasonal emission distribution in LOTOS-EUROS may require a further upward adjustment during summer. The broadly consistent pattern of change across most years in Fig. 3 supports the need to shift emissions from early spring toward the summer months to better align the model with satellite and in-situ observations.

3.2 Impact on NH_x deposition fields

Maps of the total NH_x deposition for the base and optimized simulations are shown in Fig. 4a and b. The relative-difference plots in Fig. 4c display a similar general spatial pattern as for emissions, whereby the largest increases in the modeled deposition occur in the south and east of the Netherlands, while smaller or negative differences are found in the north of the country. The largest changes in the total deposition are seen in 2020 and 2022 with differences of +10.4 % and +9.6 %, respectively. Very small relative differences are seen in 2019 and 2021 of +3.8 % and +2.7 %, respectively. The mean relative change in the NH_x deposition over the full 2018–2022 period is +6.3 %. It should be noted that a large relative difference is seen over the area of the IJsselmeer and along the Dutch coastline that is caused by the dry re-emission of NH_x , which is derived from standard maps in the LOTOS-EUROS model. However, the absolute differences (in terms of kg m⁻²) in these areas are negligible. The reduced sensitivity of deposition to the assimilation is expected, since part of the emitted NH₃ is exported offshore.

Figure 4 (a) The modeled total NH_x deposition (dry + wet) by year and the mean for the 2018–2022 period from the base simulation, (b) the same but for the optimized deposition from the LETKF assimilation run, and (c) the mean relative differences between the base and the optimized deposition fields.

3.3 Spatial distribution of NH₃ Concentrations

In this section, we evaluate the effects of the assimilation on the spatial distribution of NH₃ concentrations. Figure 5 shows the spatial distribution of NH₃ total column concentrations over the model domain from 2018 to 2022, comparing the base and optimized simulations. Panels a on the top row present the NH₃ total columns from the base model, panels b on the middle row show the updated concentration fields after integrating satellite observations, and panels c on the bottom row display the relative differences, highlighting the assimilation's impact on simulated NH₃ concentrations. Across all years, both the baseline and optimized NH₃ total column fields exhibit similar large-scale spatial patterns, with elevated NH₃ levels concentrated in regions known for intensive agricultural activity, such as northwestern Germany, and the south-eastern Netherlands. However, clear differences emerge between the prior and assimilated estimates, with the assimilation generally increasing NH₃ concentrations across most regions, as shown by the predominantly positive relative differences in panels c. The largest change in the total columns between the base and optimized simulations is seen in 2020, where a mean increase of 29.3 % is found, while the years with the smallest differences were 2019 and 2021 where an increase in the columns across the domain of roughly 3.5 % in both years was found. The mean differences over the 2018–2022 period indicate a systematic increase in mean NH₃ column concentrations within the model domain of approximately 10.4 %, with the largest increases seen over the south eastern Netherlands, suggesting that the base model likely underestimates emissions or overestimates deposition processes there. Likewise, the relative-difference panels in Fig. 5c should be interpreted in the context of the underlying base and optimized total-column fields in Fig. 5a and b, since the largest percentage changes do not always coincide with the largest absolute concentration changes. In some areas, namely in the northern and southern western parts of the domain in 2018 and 2019, slight reductions in the NH₃ total columns are seen after assimilation, which is likely due to an overestimation in the a priori emissions, and this broadly mirrors the pattern of the emissions changes from Fig. 2. For 2020 and 2022, the mean differences in NH₃ total columns are positive across the domain, whereas the corresponding emission difference maps show localized decreases, particularly in the northeastern part of the model domain. This apparent mismatch likely reflects the influence of atmospheric transport, as NH₃ emitted in one region can be advected and deposited elsewhere, so column enhancements do not necessarily coincide spatially with the emission sources. The mean of the differences for all years reveals a relatively consistent year-to-year pattern of enhancement, reinforcing the robustness of these corrections over multiple years. The systematic increases in NH₃ after assimilation suggest that satellite-derived observations are providing important constraints to correct for underestimation in the base model and emissions, particularly in regions where agricultural sources dominate. These adjustments likely reflect the assimilation compensating for structural model biases such as underestimated volatilization, overly rapid deposition, spatial misallocation in emission inventories, or the fact that the current temporal NH₃ agricultural emission profiles do not take meteorological conditions sufficiently into account (Ge et al., 2023). As discussed in Sect. 2.2, satellite retrievals can carry systematic biases arising from factors such as a priori assumptions, poor observational conditions (e.g. low thermal contrast), or cloud screening, which may influence the resulting analysis fields. By assimilating several satellite datasets, the impact of biases in individual datasets can be reduced to some extent, and this will be explored in further detail in Sect. 3.5.2. Additional support for the temporal behavior of the optimized NH₃ concentration fields is provided by the monthly MAN time series shown in Appendix Figs. B2 and B4, although these comparisons are discussed in detail later in Sect. 3.5.4.

Figure 5 (a) The mean NH₃ total column concentrations by year and for the full 2018–2022 period for the LOTOS-EUROS baseline simulation, (b) the same but for the LETKF optimized simulation, and (c) the mean relative differences between optimized and baseline simulations.

3.4 Averaging kernel sensitivity and observation density

To evaluate the impact of the number of available satellite observations on the assimilation, we examine the local observational constraint on the NH₃ emission adjustment factors using the averaging kernel sensitivity. Following Chen et al. (2023), the averaging kernel matrix is defined as

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }\mathbf{A}=\mathbf{I}-\widehat{{\mathbf{S}}^{\prime }}{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}},\end{array}$$

where I is the identity matrix, $\widehat{{\mathbf{S}}^{\prime }}$ is the posterior error covariance matrix, and S_a is the a priori error covariance matrix, both of which are directly output per time step from the LETKF. In the scalar case considered here, this reduces locally to the diagonal element $A_{ii}=\mathrm{1}-s_a/s_f$, where s_a and s_f are the local analysis and forecast error variances, respectively. The mapped values shown below therefore represent the local averaging kernel sensitivity at each grid cell and time step, expressed as percentages. Strictly speaking, the degrees of freedom for signal (DOFS) are given by the trace of A (Chen et al., 2023).

Figure 6 (a) The mean NH₃ averaging kernel sensitivity (%) in the model domain, and (b) a time series of the number of NH₃ satellite observations per month within the domain. The total number of observations per month is shown by the solid black line.

Figure 6a presents the mean spatial distribution of the averaging kernel sensitivity for each simulation year, as well as the mean over the full 2018–2022 period. Regions with broad observation coverage and high sampling density, particularly over areas with higher NH₃ concentrations, where retrieval sensitivity is generally greater and retrieval uncertainties are lower, exhibit elevated averaging kernel sensitivity values, indicating a stronger observational influence. Conversely, regions with sparse observations, such as areas with lower retrieval sensitivity, show lower values, implying a stronger dependence on the model prior. A consistent pattern of higher averaging kernel sensitivity is seen in the northwestern German border region and throughout much of the Netherlands, particularly along the western half of the country, in all years. In addition, higher mean averaging kernel sensitivity is seen in 2020, when the overlap between the CrIS and IASI instruments was greatest, as reflected in Fig. 6b. The CrIS-2 NH₃ data product begins in March 2019 and CrIS-1 observations end in mid-2021 as a result of the instrument being decommissioned, while the IASI-C NH₃ data product starts in September 2019 and the IASI-A product ended in October 2021 due to the MetOp-A platform reaching the end of its operational lifetime.

To assess the relationship between averaging kernel sensitivity and observation density, a time series of the total number of observations in the model domain is shown in Fig. 6b. The higher mean averaging kernel sensitivity in 2020 coincides with the period of greatest overlap in CrIS and IASI availability, suggesting that increased observation coverage contributed to stronger observational constraint in that year. This impact is also consistent with Figs. 2 and 5, where the largest adjustments in emissions and total column concentrations relative to the original simulation were observed in 2020. It is possible that in years with more limited observations, the post-analysis emission changes are underestimated because less observational information is available to constrain the system. However, as discussed earlier, 2020 was also an exceptional year with very warm summer months and particularly high NH₃ emissions, which likely also contributed to the higher averaging kernel sensitivity. The relationship is therefore not strictly linear, since the averaging kernel sensitivity depends not only on observation count, but also on the prior and posterior error covariance structure and on the relative uncertainty of the satellite retrievals. This means that the assimilation can be limited both by the density of available retrievals and by their accuracy: high observational coverage with large uncertainties provides limited benefit, whereas accurate retrievals at low coverage can still meaningfully constrain the system, but are less able to capture transient emission events or resolve daily and diurnal cycles. In practice, both factors act together, and the results emphasize the importance of maintaining dense, accurate satellite observations, ideally from multiple instruments, to more effectively constrain NH₃ emissions and capture interannual variability.

3.5 Comparisons with ground-based observations

Although the assimilation produces consistent adjustments to the modeled emissions, concentrations, and to a lesser extent the deposition fields across the simulation period, independent evaluation against ground-based observations is needed to assess the extent to which these optimized fields represent an improvement over the original simulation.

3.5.1 LML hourly surface concentrations

To accurately assess the impact of the assimilation of the satellite observations, and to evaluate whether the optimized concentration fields represent more realistic estimates of the true spatio-temporal distribution of NH₃, it is important to compare the results against an independent dataset. The surface concentrations from the base and optimized model runs were compared against ground-based surface observations from six sites in the LML network that were described in Sect. 2.3.1. These sites were selected because they provide hourly NH₃ surface concentration measurements and have significant time-series (i.e. >1 year) of data. As the primary focus of this study is NH₃, detailed discussion of the corresponding NO₂ assimilation results is beyond the scope of the present study. These results will be presented separately in a forthcoming paper to allow for a more comprehensive treatment.

A correlation plot of the monthly temporal comparison is shown in Fig. 7, and a corresponding plot of the spatial means (i.e. all sites averaged for a given month) is provided in Fig. 8. These two summary views are shown separately to distinguish temporal agreement at individual sites from agreement in the spatial pattern across the Dutch monitoring network. They provide compact statistical summaries, while complementary temporal and spatial context is given by the diurnal-cycle analysis and the mapped concentration fields shown elsewhere in the manuscript. All linear regressions were performed using an ordinary least squares fitting approach. Uncertainty estimates are reported for the statistical quantities shown in the model–observation scatter-plot comparisons. Unless otherwise stated, uncertainties correspond to one standard error. For the mean bias, μ, the uncertainty was estimated as $s_d/\sqrt{N}$, where s_d is the sample standard deviation of the model–observation differences and N is the number of paired data points. For the standard deviation of the differences, σ, the standard error was approximated as $\mathit{\sigma }/\sqrt{\mathrm{2}(N-\mathrm{1})}$, assuming normally distributed differences. For the Pearson correlation coefficient, R, uncertainties were estimated using the Fisher transformation, with standard error $\mathrm{1}/\sqrt{N-\mathrm{3}}$ in Fisher-z space and then transformed back to correlation space. Uncertainties in linear regression slopes were estimated from the ordinary least-squares covariance matrix. These uncertainty estimates assume independent paired samples and should therefore be interpreted as approximate, particularly for comparisons involving repeated monthly values across sites. For changes between the base and optimized simulations, statistical robustness was assessed using paired bootstrap resampling of the matched observation, base-model, and optimized-model triplets; changes were considered statistically significant when the 95 % bootstrap confidence interval did not include zero. The main surface NH₃ comparison statistics are summarized in Table 2, while the full surface-concentration uncertainty estimates and bootstrap confidence intervals are provided in Appendix Tables C1 and C2. Corresponding uncertainty estimates and bootstrap confidence intervals for the wet deposition comparisons are provided in Appendix Tables C3 and C4. Uncertainty estimates for the 2020 satellite-subset sensitivity experiments are provided in Appendix Table C5.

Figure 7 Scatter plot of monthly temporal means of (left) LOTOS-EUROS base NH₃, and (right) LOTOS-EUROS LETKF optimized NH₃ vs. LML observed NH₃ surface concentrations for the period of January 2018–December 2022. Uncertainty estimates for the reported statistics and paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Tables C1 and C2.

Table 2 Summary of model performance statistics for NH₃ surface concentration comparisons against LML and MAN observations for the 2018–2022 period. Reported values correspond to monthly spatial-mean comparisons. Values for the slope, mean bias (μ), and spread of model–observation differences (σ) are reported as estimate ± one standard error. Confidence intervals for R correspond to the 95 % interval obtained using the Fisher transformation.

Not all LML sites provided data for all months, leading to a total of 346 data points compared with the expected 360 in Fig. 7. For the LML temporal comparison, the assimilation produces a statistically robust reduction in mean bias and a significant improvement in the regression slope, while the already high temporal correlation remains statistically similar before and after assimilation. The bias decreases from $-\mathrm{1.97}\pm \mathrm{0.21}$ to $-\mathrm{0.97}\pm \mathrm{0.21}$ µg m⁻³, and paired bootstrap resampling indicates that this change is significant at the 95 % level ($\mathrm{\Delta }\mathit{\mu }=+\mathrm{0.99}$ µg m⁻³, 95 % CI: +0.74–+1.25 µg m⁻³). The regression slope increases from 0.790±0.028 to 0.909±0.031, also representing a statistically significant improvement ($\mathrm{\Delta }\text{slope}=+\mathrm{0.118}$, 95 % CI: +0.062–+0.173). In contrast, the change in correlation from R=0.836 [0.801, 0.865] to R=0.846 [0.813, 0.874] is not statistically significant ($\mathrm{\Delta }R=+\mathrm{0.015}$, 95 % CI: −0.012–+0.040). For the monthly spatial means shown in Fig. 8, the assimilation again produces a statistically significant reduction in mean bias, from $-\mathrm{2.04}\pm \mathrm{0.36}$ to $-\mathrm{1.05}\pm \mathrm{0.31}$ µg m⁻³ ($\mathrm{\Delta }\mathit{\mu }=+\mathrm{0.99}$ µg m⁻³, 95 % CI: +0.54–+1.44 µg m⁻³). The Pearson correlation and regression slope also increase, from R=0.765 [0.635, 0.853] to R=0.837 [0.740, 0.900] and from 0.787±0.087 to 0.895±0.077, respectively. However, paired bootstrap confidence intervals for these changes include zero ($\mathrm{\Delta }R=+\mathrm{0.071}$, 95 % CI: −0.002–+0.146; $\mathrm{\Delta }\text{slope}=+\mathrm{0.110}$, 95 % CI: −0.028–+0.236), so these increases should be interpreted as suggestive rather than statistically significant at the 95 % level. The LML comparison therefore indicates that the most robust improvement is the reduction of systematic underestimation, with additional but weaker evidence for improved spatial agreement.

Figure 8 Scatter plot of monthly spatial means of (left) LOTOS-EUROS base NH₃, and (right) LOTOS-EUROS LETKF optimized NH₃ vs. LML observed NH₃ surface concentrations for the period of January 2018–December 2022. Each data-point represents the mean calculated across all LML sites for a given month, and are colored corresponding to the month while the marker style indicates the year. Uncertainty estimates for the reported statistics and paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Tables C1 and C2.

A Taylor diagram illustrating the measurement-model comparisons for the base and optimized simulations at each individual LML site is shown in Fig. 9. In general, the correlations between the measured and modeled surface NH₃ concentrations are higher for the optimized run at all LML sites except Wekerom-Riemterdijk, with the clearest increases seen at Zegveld Oude-Meije and Valthermond-Noorderdiep. The correlation shows a small improvement at Vredepeel-Vredeweg and this was accompanied by a substantial reduction in the bias from −3.10 to 0.20 µg m⁻³, but a large increase in the relative standard deviation is also seen that is driven by an apparent larger number of short-term enhancement events being captured in the modeled time-series in the optimized simulation (shown in Appendix Fig. A1). The mean over all sites displays a small improvement in the correlation, but also a small increase in the standard deviation (also shown earlier in Fig. 7).

Figure 9 Taylor diagram displaying the Pearson correlation coefficients and the relative standard deviations of the comparisons between the model and each LML site for the 2018–2022 period for (black) the baseline simulation and (red) the optimized simulation.

Figure 10 Diurnal distributions of surface NH₃ 2018–2022 at each LML site from the observations (dark grey), the base model simulation (red), and the optimized simulation (blue). The box-and-whisker representation shows the distribution of hourly NH₃ concentrations at each hour of the day for the three datasets. On each panel, the difference between the median diurnal cycles (Δ_median) and the Pearson correlation coefficient (R) are provided for the base run in red and the optimized run in blue.

Since the LML data are provided at an hourly frequency, the impact of assimilation on the diurnal cycles of NH₃ in the model can also be investigated at each site. The mean diurnal cycles from the observations, the base model simulation, and the optimized simulation calculated over the 2018–2022 period are shown in Fig. 10. The box-and-whisker representation is used to show not only the central tendency of the diurnal cycle, but also the spread of the hourly concentration distributions at each site. In most cases, the mean differences in the diurnal cycles showed improvement relative to the observations in the optimized run in comparison to the base simulation even though only morning and afternoon satellite overpasses were used. This is partly because the effect of the assimilation persists between overpass times through the forecast step and temporal persistence of the emission adjustments, allowing the updated state to influence concentrations beyond the observation times themselves. From Fig. 10 it can be seen that in the base simulation, the diurnal cycles are in many cases largely underestimated relative to the surface observations, however after the assimilation of the satellite observations a much closer agreement between the measurements and the model is found. At Wekerom-Riemterdijk and Wieringerwerf-Medemblikkerweg the mean differences increased from −2.03 to −2.70 and −2.83 to −2.92 µg m⁻³, respectively. At all other LML sites a decrease in the mean difference of the diurnal cycles was observed in the optimized run relative to the base run. The correlations in the mean diurnal cycles remain largely the same, with the exception of Wieringerwerf-Medemblikkerweg where a moderate improvement from R=0.72 to R=0.77 is found. However, the diurnal cycle at Wieringerwerf-Medemblikkerweg appears to be substantially underestimated in both the base model and the optimized simulation. These underestimations are likely related to local-scale NH₃ enhancements from the nearby farmlands in the Wieringermeer polder being poorly captured by the model, coupled with complex local meteorological conditions due to the proximity of the site to the IJselmeer and the North Sea which is roughly 15–20 km away.

Although the assimilation of satellite observations improves the representation of NH₃ diurnal cycles relative to measurements at several LML sites, particularly with respect to systematic biases, it does not fully reproduce the observed variability. A fundamental limitation of Kalman filter–based approaches is that, although emissions are adjusted, there can be a temporal lag associated with the transport of concentrations from source regions to the measurement locations, which is strongly modulated by local meteorological conditions. Furthermore, when the underlying diurnal cycle in the model is misrepresented, assimilation at only two main satellite overpass times may be insufficient to correct these deficiencies. An additional challenge is that NH₃ emissions originate from diverse sources such as livestock housing, manure application, and industry, each with distinct diurnal patterns, whereas the present assimilation setup applies a single total emission adjustment per species and grid cell, without distinguishing among source types that may have different underlying diurnal emission profiles. Overall, at sites located in regions of intensive agricultural activity such as Valthermond, Vredepeel, and Zegveld where midday enhancements are more likely to be captured by satellite observations, the influence of the assimilation on the diurnal cycles is more pronounced.

While these comparisons demonstrate clear improvements after assimilation, several further limitations should be acknowledged. First, the six LML sites provide only sparse coverage of the Netherlands and may not fully represent the strong spatial heterogeneity of NH₃ concentrations across the entire country, particularly in regions with intensive agricultural emissions. Second, the model–observation comparison involves a scale mismatch: the LML instruments measure point concentrations, while the LOTOS-EUROS model output represents grid-cell averages at 7 km×7 km resolution, which can potentially introduce representativity errors. In addition, the ground-based measurements themselves are not perfect and subject to calibration uncertainties and potential interferences. Lastly, because the statistics are based on monthly averages, the analysis does not explicitly resolve short-term (daily to synoptic) variability. Consequently, part of the observed agreement may reflect the influence of meteorological variability in addition to the direct effects of the assimilation. Taken together, the LML comparisons suggest that the assimilation generally improves the agreement of the model with independent ground-based observations, particularly through a robust reduction in systematic underestimation, with weaker but positive evidence for improved spatial agreement. At the same time, the analysis also underlines the limitations of both the observational dataset and the comparison methodology, which should be taken into account when interpreting the results.

3.5.2 Impact of satellite selection on LML surface comparisons

To investigate the impact of the co-assimilation and the choice of satellites on the final optimized model state, assimilation runs were repeated for the year of 2020 using subsets of the satellite data products. These runs were: IASI only, CrIS only, and CrIS and IASI (without TROPOMI NO₂), and the comparisons against the hourly LML NH₃ observations were then repeated.

The temporal means for each of the four assimilation runs vs. the LML observations are shown in Fig. 11. From Fig. 11, it can be seen that relative to the other assimilation runs, the assimilation performed using IASI, CrIS and TROPOMI yields the lowest overall bias (0.1 µg m⁻³) relative to the LML observations while maintaining a relatively strong correlation of R=0.88. The simulation performed with IASI and CrIS without the assimilation of TROPOMI NO₂ yields a higher correlation of R=0.90, but displays a higher bias of 0.7 µg m⁻³, a higher standard deviation, and a poorer regression slope. The simulations using only CrIS and only IASI display relatively higher mean biases of 1.4 and −1.1 µg m⁻³, respectively. For comparison purposes, the temporal mean comparison for the base (unoptimized) simulation for 2020 is shown in the Appendix in panel a of Fig. A2. Each of the simulations performed with assimilation show an improvement over the unoptimized base model run, which displays a correlation of R=0.82 and a mean bias of −1.9 µg m⁻³ relative to the LML observations.

Figure 11 Scatter plot of monthly temporal means of assimilation runs using (a) IASI and CrIS NH₃ and TROPOMI NO₂ observations, (b) IASI and CrIS NH₃ observations, (c) CrIS NH₃ observations only, and (d) IASI NH₃ observations only vs. LML observed NH₃ surface concentrations for the year of 2020. Uncertainty estimates for the reported statistics are provided in Appendix Table C5.

A similar plot of the spatial means for each of the assimilation runs vs. the LML surface observations is provided in Fig. 12. It can be seen that again from Fig. 12a, the assimilation performed using IASI, CrIS, and TROPOMI leads to the lowest overall bias (0.1 µg m⁻³) of the four simulations and a relatively strong correlation of R=0.90. In comparison, the assimilation run using IASI and CrIS observations but not TROPOMI has a slightly higher correlation (R=0.93), but a higher bias of 0.7 µg m⁻³ and a poorer slope of regression (1.163 vs. 0.977 from the IASI, CrIS and TROPOMI assimilation). The assimilation runs performed using CrIS only and IASI only (shown in panels c and d of Fig. 12, respectively) have the same or poorer correlations than the IASI, CrIS and TROPOMI run, and relatively larger biases of 1.4 and −1.1 µg m⁻³ for the CrIS-only and IASI-only runs, respectively. As was the case for the temporal mean comparisons, each of the optimized runs shows better agreement with the observations than the unoptimized simulation (panel b of Fig. A2), which had a weaker correlation of R=0.66 and a mean bias of −1.9 µg m⁻³ for 2020. Uncertainty estimates for the statistics reported in Figs. 11 and 12, together with the corresponding 2020 base-model comparison in Appendix Fig. A2, are provided in Appendix Table C5. Because these experiments are intended as a sensitivity analysis rather than a formal pairwise model-selection test, we report uncertainty estimates for each run but do not assess all pairwise differences between satellite configurations.

Figure 12 Scatter plot of monthly spatial means of assimilation runs using (a) IASI and CrIS NH₃ and TROPOMI NO₂ observations, (b) IASI and CrIS NH₃ observations, (c) CrIS NH₃ observations only, and (d) IASI NH₃ observations only vs. LML observed NH₃ surface concentrations for the year of 2020. Uncertainty estimates for the reported statistics are provided in Appendix Table C5.

These comparisons highlight that the choice of satellite datasets used for assimilation is impactful, and that using just IASI or CrIS for the optimization may lead to greater biases than using a combined assimilation approach. Although assimilating a combination of CrIS and IASI observations leads to an improvement over just a single subset of the satellite data, the inclusion of TROPOMI NO₂ in the optimization aids in reducing the overall biases while still maintaining a strong correlation, particularly in comparison to the un-optimized LOTOS-EUROS simulation. The added value of assimilating NO₂ is that it helps to constrain the availability of HNO₃ through the NO₂–OH oxidation pathway, which in turn directly regulates the partitioning of NH₃ into ammonium nitrate, an important sink of NH₃ (Finlayson-Pitts and Pitts, 2000). In this way, the NO₂ observations provide an indirect but important constraint on the fate of NH₃, thereby improving the consistency of the coupled reactive nitrogen system in the model. Overall, the subset experiments indicate that the multi-satellite configuration including IASI, CrIS, and TROPOMI provides the best balance between low bias, strong correlation, and physically consistent reactive-nitrogen constraints, although the sensitivity experiments are not treated here as a formal pairwise model-selection test.

3.5.3 LML NH_x wet deposition comparisons

To evaluate the impact of the assimilation on the deposition of NH_x within the Netherlands, comparisons were performed against dissolved ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations measured in precipitation. These measurements are made at irregular intervals at several LML sites, and for the comparisons, these observations were paired with the modeled deposition interpolated to a daily frequency and converted to wet deposition fluxes (in $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) using the measured precipitation amounts.

Figure 13 Scatter plot of monthly temporal means of (left) base LOTOS-EUROS ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition flux, and (right) LOTOS-EUROS LETKF optimized ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition flux vs. LML surface observations for the period of January 2018–December 2022. Uncertainty estimates for the reported statistics and paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Tables C3 and C4.

Figure 13 shows the comparison of the measured and modeled monthly temporal means of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition flux before and after assimilation. The baseline simulation displays a moderate correlation with the observations, increasing from R=0.510 [0.410, 0.598] to R=0.560 [0.466, 0.641] after assimilation. Paired bootstrap resampling indicates that this increase is statistically significant, although modest in magnitude ($\mathrm{\Delta }R=+\mathrm{0.044}$, 95 % CI: +0.019–+0.071). The mean bias is reduced from $-\mathrm{1.00}\pm \mathrm{0.28}$ to $-\mathrm{0.80}\pm \mathrm{0.25}$ $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, and this reduction is also statistically significant ($\mathrm{\Delta }\mathit{\mu }=+\mathrm{0.21}$ $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, 95 % CI: +0.07–+0.34 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). The spread of the model–observation differences decreases significantly from 4.40±0.20 to 3.90±0.18 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ ($\mathrm{\Delta }\mathit{\sigma }=-\mathrm{0.51}$ $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, 95 % CI: −0.81–−0.21 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). The regression slope remains statistically similar before and after assimilation ($\mathrm{\Delta }\text{slope}=-\mathrm{0.009}$, 95 % CI: −0.068–+0.046). The scatter plot also shows a particularly large overestimation in the modeled deposition flux for several LML sites in the spring months of 2018 and 2020, which is reduced after assimilation, while there is a general underestimation during much of the remainder of the year. A similar pattern of springtime overestimation and late-year underestimation was found in the wet deposition comparisons in the earlier LOTOS-EUROS LETKF assimilation study by Van Der Graaf et al. (2022), who attributed this behavior to a potential overestimation of springtime NH₃ emissions and underestimation later in the year.

A scatterplot of the monthly spatial means is shown in Fig. 14. For this comparison, the correlation increases from R=0.578 [0.380, 0.725] to R=0.684 [0.520, 0.799], and paired bootstrap resampling indicates that this increase is significant at the 95 % level ($\mathrm{\Delta }R=+\mathrm{0.104}$, 95 % CI: +0.057–+0.166). The spread of the model–observation differences is also significantly reduced, from 3.04±0.28 to 2.38±0.22 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ ($\mathrm{\Delta }\mathit{\sigma }=-\mathrm{0.62}$ $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, 95 % CI: −1.08–−0.20 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$). The mean bias is reduced from $-\mathrm{1.15}\pm \mathrm{0.39}$ to $-\mathrm{0.91}\pm \mathrm{0.31}$ $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, although this change is not statistically significant at the 95 % level.

Figure 14 Scatter plot of monthly spatial means of (left) base LOTOS-EUROS ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition flux, and (right) LOTOS-EUROS LETKF optimized ${\mathrm{NH}}_{\mathrm{4}}^{+}$ wet deposition flux vs. LML surface observations for the period of January 2018–December 2022. Each data-point represents the mean calculated across all LML sites for a given month. Uncertainty estimates for the reported statistics and paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Tables C3 and C4.

Overall, the wet deposition evaluation provides a complementary downstream consistency check. The optimized simulation shows statistically significant improvements in the correlation, mean bias, and spread of the monthly temporal-mean comparison, as well as in the correlation and spread of the monthly spatial-mean comparison. The monthly spatial mean bias is also reduced, although this change is not statistically significant at the 95 % level. These results indicate that the assimilation improves several aspects of the modeled deposition field, but that wet deposition remains a less direct and less sensitive validation constraint than the surface NH₃ concentration networks.

3.5.4 MAN monthly NH₃ comparisons

The baseline and optimized simulations were also evaluated against monthly MAN NH₃ surface observations from a total of 315 sites (309 standard sites, and 6 MAN instruments located at LML sites) across the Netherlands for 2018–2022. Figure 15 shows the mean observed and modeled NH₃ concentrations, along with their differences. The baseline simulation exhibits a clear spatial bias, underestimation along the western coast and overestimation in the east, while the optimized run amplifies this pattern, consistent with the emission and concentration increases in the southern and eastern Netherlands. No comparable bias pattern was observed for the hourly LML sites or the CrIS and IASI total column comparisons.

Figure 15 Spatial maps of the mean NH₃ surface concentrations from (a) MAN observations, (b) the base simulation, and (c) the optimized simulation calculated over the period of 2018–2022. The differences (in µg m⁻³) are shown in the second row for (d) the optimized run vs. the base run, (e) the base run vs. the observations, and (f) the optimized run vs. the observations, and the relative differences (in %) are shown for the same comparisons in panels (g) to (i).

Figure 16 compares monthly spatial means across all MAN sites, and for comparison purposes, the LML and MAN statistics are summarized together in Table 2. The assimilation substantially improves the spatial correlation, from R=0.791 [0.672, 0.870] to R=0.892 [0.824, 0.934]. Paired bootstrap resampling confirms that this increase is statistically significant ($\mathrm{\Delta }R=+\mathrm{0.100}$, 95 % CI: +0.036–+0.178). However, the improved spatial correlation is accompanied by a statistically significant increase in the positive mean bias, from $+\mathrm{1.22}\pm \mathrm{0.37}$ to $+\mathrm{2.43}\pm \mathrm{0.34}$ µg m⁻³ ($\mathrm{\Delta }\mathit{\mu }=+\mathrm{1.22}$ µg m⁻³, 95 % CI: +0.71–+1.73 µg m⁻³), and by an increase in the regression slope from 1.385±0.141 to 1.651±0.110 ($\mathrm{\Delta }\text{slope}=+\mathrm{0.270}$, 95 % CI: +0.033–+0.515). The change in the spread of the model–observation differences is not statistically significant ($\mathrm{\Delta }\mathit{\sigma }=-\mathrm{0.17}$ µg m⁻³, 95 % CI: −0.82–+0.49 µg m⁻³).

Figure 16 Scatter plot of monthly spatial means of (left) base LOTOS-EUROS NH₃, and (right) LOTOS-EUROS LETKF optimized NH₃ vs. MAN observed NH₃ surface concentrations for the period of January 2018–December 2022. Each data-point represents the mean calculated across all MAN sites for a given month, and are colored corresponding to the month while the marker style indicates the year. Uncertainty estimates for the reported statistics are provided in Table 2 and Appendix Table C1; paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Table C2.

These results indicate that the LETKF-optimized simulation better captures the large-scale spatial variability sampled by the MAN network, but that this improved spatial coherence does not correspond to an overall reduction in model error for the full MAN dataset. Instead, the optimized simulation increases the positive offset relative to MAN observations. This contrasts with the LML surface NH₃ comparison, where the assimilation significantly reduced the negative bias and brought the regression slope closer to unity. The different behavior of the full MAN network is consistent with the representativeness interpretation discussed above: many MAN sites are located in Natura 2000 areas, often in sheltered or forested environments where local NH₃ concentrations are lower than the surrounding landscape-scale concentrations represented by the 7 km×7 km model grid and by the satellite footprints used in the assimilation. Supporting visual comparisons using individual-site MAN monthly values and national monthly mean time series (Appendix Figs. B1 and B2) show a similar qualitative pattern to the full MAN-network comparison, with the optimized run better reproducing broad interannual variability but maintaining a high bias, particularly during the summer months of 2020 and 2022, the years with the strongest emission adjustments (see Fig. 2).

To further assess the spatial representativity of the MAN network and to enable a more direct comparison with the LML results, we also evaluated the six MAN calibration sensors co-located with LML sites against the base and optimized LETKF simulations. These six sensors are not part of the standard MAN dataset, but are used for the monthly calibration of the remaining 309 MAN sites, as described by Noordijk et al. (2020). A spatial-mean scatterplot for these six sensors is provided in Appendix Fig. B3, and the associated statistics are summarized in Table 2. In contrast to the full MAN network, the co-located MAN sensors show behavior consistent with the LML analysis. The mean bias is significantly reduced from $-\mathrm{2.42}\pm \mathrm{0.46}$ to $-\mathrm{0.99}\pm \mathrm{0.41}$ µg m⁻³ after assimilation ($\mathrm{\Delta }\mathit{\mu }=+\mathrm{1.42}$ µg m⁻³, 95 % CI: +0.85–+2.01 µg m⁻³). The spatial correlation also improves significantly from R=0.662 [0.490, 0.784] to R=0.761 [0.629, 0.851] ($\mathrm{\Delta }R=+\mathrm{0.098}$, 95 % CI: +0.011–+0.191), and the regression slope increases from 0.666±0.099 to 0.831±0.093 ($\mathrm{\Delta }\text{slope}=+\mathrm{0.167}$, 95 % CI: +0.012–+0.325). The spread of the model–observation differences decreases from 3.57±0.33 to 3.15±0.29 µg m⁻³, although this change is not statistically significant. The much better agreement obtained for the co-located MAN calibration sensors indicates that the broader positive bias in the full MAN network is dominated primarily by representativeness differences rather than by a uniform model overestimation across all MAN locations. In other words, when the MAN comparison is restricted to more open, regionally representative sites similar to the LML locations, the optimized simulation shows statistically robust improvements in bias, slope, and spatial correlation. Part of this contrast may stem from measurement uncertainties in the passive Gradko samplers, but representativeness effects are expected to dominate. The LML stations and the satellite footprints sample air masses that are more representative of the broader agricultural landscape, whereas most MAN sites are located within Natura 2000 areas, frequently in or near forested environments where canopy shielding and strong local heterogeneity depress NH₃ concentrations relative to the surrounding terrain. These sheltered conditions are highly localized and not captured by the coarse spatial scales of the regional model (7 km×7 km) or the 12–14 km nadir footprint of the satellite retrievals. As a result, the MAN measurements at these locations tend to reflect sub-grid processes rather than the landscape-average conditions that the model and satellite products are designed to represent. Remaining differences may also relate to how canopy deposition or near-surface mixing is expressed at these fine spatial scales, rather than broad model biases.

Overall, the MAN comparisons reaffirm that the LETKF enhances the spatial coherence of the NH₃ fields and performs consistently when evaluated against regionally representative measurements such as the LML sites. The discrepancies observed at many MAN locations primarily reflect representativeness limitations inherent in comparing coarse-scale model and satellite products with low-temporal-resolution passive samplers situated in highly heterogeneous environments. Reducing these mismatches will likely require finer-scale process representation and more explicit treatment of land-cover–dependent effects in future assimilation studies.

4 Conclusions

This study demonstrates that assimilating NH₃ and NO₂ satellite observations into the LOTOS-EUROS model using the LETKF framework substantially improves the representation of reactive nitrogen dynamics over the Netherlands. By co-assimilating IASI, CrIS, and TROPOMI retrievals over 2018–2022, the system produced optimized emissions, deposition, and concentration fields. The optimized emission fields showed consistent spatial structures across years, with persistent increases in the southern and eastern Netherlands, and exhibited notable temporal shifts such as a reduced springtime emission peak and enhanced summertime emissions. These patterns are broadly consistent with recent observational and modeling studies evaluating NH₃ emissions over the region (Ge et al., 2020; Van Der Graaf et al., 2022; Ding et al., 2024).

Validation against independent LML surface observations at six sites showed a statistically robust reduction in mean bias, from approximately −2.0 to −1.0 µg m⁻³, while changes in correlation were positive but more modest. Assimilation also improved the representation of diurnal cycles, reducing systematic underestimation and bringing simulated amplitudes closer to observations, although some residual site-specific variability remained. Comparisons with dissolved ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in precipitation provided a complementary downstream evaluation, showing statistically significant improvements in temporal correlation, temporal bias, temporal spread, and in the spatial correlation and spread of wet deposition differences, while the monthly spatial mean-bias change was more modest and not statistically significant. Sensitivity tests for 2020 indicated that the multi-satellite configuration performed best overall, with joint assimilation of IASI and CrIS NH₃ and TROPOMI NO₂ producing the lowest biases relative to LML.

Comparisons with monthly MAN observations showed a statistically significant improvement in spatial correlation after assimilation, from R=0.79 to R=0.89, but also a statistically significant increase in positive bias, from +1.2 to +2.4 µg m⁻³, for the full MAN network. These biases primarily reflect representativeness differences: many MAN sites are located within Natura 2000 areas, often in forested or sheltered environments where canopy effects and fine-scale heterogeneity suppress local NH₃ levels relative to the surrounding agricultural landscape. Such localized variability is not resolved at the model grid scale or by satellite footprints, leading to systematic overestimation when landscape-scale fields are evaluated against sheltered MAN locations.

In contrast, the MAN calibration sensors co-located at six LML stations showed statistically robust improvements in bias, slope, and spatial correlation. At these sites, the pre-assimilation mean bias of −2.4 µg m⁻³ improved to −1.0 µg m⁻³, the regression slope increased from 0.67 to 0.83, and the spatial correlation increased from R=0.66 to R=0.76. These results support the interpretation that the broader MAN bias is primarily driven by representativeness differences. Together, the LML and MAN comparisons highlight both the strengthened large-scale spatial and temporal structure of NH₃ after assimilation and the ongoing challenges of reconciling coarse-resolution model and satellite fields with passive samplers located in highly heterogeneous environments. Because the present LETKF setup optimizes total NH₃ emissions rather than sector-resolved source contributions, the spatial emission adjustments shown here cannot be attributed robustly to individual source types. Further progress may be enabled by adopting a label-based Kalman filtering approach. The labeling functionality introduced in LOTOS-EUROS v2.3 could be extended to the LETKF, allowing sector-specific emission optimization and supporting finer-scale improvements.

The results also underscore the importance of satellite observation density, particularly for NH₃. The year 2020, characterized by the highest availability of CrIS and IASI retrievals, showed the largest emission and concentration adjustments. The recent launch of MTG-IRS on the MTG-S1 platform (1 July 2025) presents a major opportunity: its expected half-hourly 4 km×4 km NH₃ retrieval capability could provide unprecedented spatial and temporal coverage over Europe, enabling improved constraints on diurnal variability, better characterization of emission events, and more rigorous evaluation of emission inventories. A dedicated NH₃ assimilation study using MTG-IRS is recommended once data become available.

In addition to these advancements, several refinements warrant further investigation. These include improving the representativeness between model and surface observations, enhancing the treatment of diurnal variability across source sectors, and conducting a full system simulation experiment (SSE) to rigorously evaluate LETKF performance.

Overall, these findings highlight the essential role of satellite constraints in advancing chemical transport modeling and nitrogen budget estimation. Co-assimilating complementary satellite instruments provides a pathway toward more accurate and internally consistent representations of reactive nitrogen, strengthening our ability to constrain emissions, evaluate inventories, and assess deposition at regional to national scales. The demonstrated improvements offer a strong foundation for exploiting emerging missions such as MTG-IRS, enabling continued progress toward capturing fine-scale processes and further improving assimilation performance in future studies.

Appendix A: Additional figures for LML network comparisons

A1 Vredepeel-Vredeweg hourly NH₃ timeseries

Figure A1 Timeseries of hourly NH₃ surface concentration values at the Vredepeel-Vredeweg LML site from (black) the observations, (red) the base LOTOS-EUROS simulation, and (blue) the LETKF-optimized simulation.

A2 LML comparisons for 2020 with base model simulation

Figure A2 Scatter plot of (a) monthly temporal means of LOTOS-EUROS simulated NH₃ vs. LML observed NH₃ surface concentrations, and (b) monthly spatial means of LOTOS-EUROS simulated NH₃ vs. LML observed NH₃ surface concentrations for the period of January–December 2020. Uncertainty estimates for the reported statistics are provided in Appendix Table C5.

Appendix B: Additional figures for MAN network comparisons

B1 Monthly mean scatter plot pre- and post-assimilation

Figure B1 Scatter plot comparing monthly means of observed MAN NH₃ mass concentrations with (left) base LOTOS-EUROS surface NH₃ mass concentration, and (right) LOTOS-EUROS LETKF optimized NH₃ mass concentration. Each data-point represents a mean at a single MAN site. Uncertainty estimates for the reported statistics and paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Tables C1 and C2.

B2 Monthly mean timeseries

Figure B2 Time-series of monthly mean NH₃ surface mass concentrations and differences calculated across all MAN sites during 2018–2022. The shaded regions indicate the standard deviations of the monthly means. The mean differences shown in the monthly mean time series correspond to the same monthly spatial-mean comparisons summarized in Table 2 and Appendix Tables C1 and C2.

B3 MAN @ LML sites: spatial correlations

Figure B3 Scatter plot of monthly spatial means of (left) base LOTOS-EUROS NH₃, and (right) LOTOS-EUROS LETKF optimized NH₃ vs. MAN observed NH₃ surface concentrations at the locations of LML sites for the period of January 2018–December 2022. Each data-point represents the mean calculated across the 6 MAN calibration sensors located at the LML sites for a given month, and are colored corresponding to the month while the marker style indicates the year. Uncertainty estimates for the reported statistics are provided in Table 2 and Appendix Table C1; paired bootstrap confidence intervals for optimized-minus-base changes are provided in Appendix Table C2.

B4 MAN @ LML sites: monthly mean time-series

Figure B4 Time-series of monthly mean NH₃ surface mass concentrations and differences calculated across the 6 MAN calibration sensors located at the LML sites during 2018–2022. The shaded regions indicate the standard deviations of the monthly means. The mean differences shown in the monthly mean time series correspond to the same monthly spatial-mean comparisons summarized in Table 2 and Appendix Tables C1 and C2.

Appendix C: Detailed uncertainty estimates for model–observation statistics

C1 Uncertainty estimates for surface NH₃ comparisons

Table C1 Full uncertainty estimates for model–observation comparison statistics. Values for the slope, mean bias (μ), and spread of model–observation differences (σ) are reported as estimate ± one standard error. Confidence intervals for R correspond to the 95 % interval obtained using the Fisher transformation. Units for μ and σ are µg m⁻³.

C2 Bootstrap assessment of optimized-minus-base changes for surface NH₃ comparisons

Table C2 Paired bootstrap estimates of optimized-minus-base changes in model–observation statistics. For each case, bootstrapping was performed with 10 000 bootstrap resamples. Confidence intervals represent 95 % bootstrap intervals. Statistically significant changes are those for which the confidence interval does not include zero. Units for Δμ and Δσ are µg m⁻³.

C3 Uncertainty estimates for wet deposition comparisons

Table C3 Uncertainty estimates for wet deposition model–observation comparison statistics. Values for the slope, mean bias (μ), and spread of model–observation differences (σ) are reported as estimate ± one standard error. Confidence intervals for R correspond to the 95 % interval obtained using the Fisher transformation. Units for μ and σ are $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$.

C4 Bootstrap assessment of optimized-minus-base changes for wet deposition comparisons

Table C4 Paired bootstrap estimates of optimized-minus-base changes in wet deposition model–observation statistics. For each case, bootstrapping was performed with 10 000 bootstrap resamples. Confidence intervals represent 95 % bootstrap intervals. Statistically significant changes are those for which the confidence interval does not include zero. Units for Δμ and Δσ are $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$.

C5 Uncertainty estimates for satellite-subset sensitivity experiments

Table C5 Uncertainty estimates for the 2020 satellite-subset sensitivity experiments evaluated against LML NH₃ surface observations. Values for the slope, mean bias (μ), and spread of model–observation differences (σ) are reported as estimate ± one standard error. Confidence intervals for R correspond to the 95 % interval obtained using the Fisher transformation. Units for μ and σ are µg m⁻³.