NAQPMS-PDAF v2.0: a novel hybrid nonlinear data assimilation system for improved simulation of PM_2.5 chemical components

Identifying PM_2.5 chemical components is crucial for formulating emission strategies, estimating radiative forcing, and assessing human health effects. However, accurately describing spatiotemporal variations in PM_2.5 chemical components remains a challenge. In our earlier work, we developed an aerosol extinction coefficient data assimilation (DA) system (Nested Air Quality Prediction Model System with the Parallel Data Assimilation Framework (NAQPMS-PDAF) v1.0) that was suboptimal for chemical components. This paper introduces a novel hybrid nonlinear chemical DA system (NAQPMS-PDAF v2.0) to accurately interpret key chemical components (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, OC, and EC). NAQPMS-PDAF v2.0 improves upon v1.0 by effectively handling and balancing stability and nonlinearity in chemical DA, which is achieved by incorporating the non-Gaussian distribution ensemble perturbation and hybrid localized Kalman–nonlinear ensemble transform filter with an adaptive forgetting factor for the first time. The dependence tests demonstrate that NAQPMS-PDAF v2.0 provides excellent DA results with a minimal ensemble size of 10, surpassing previous reports and v1.0. A 1-month DA experiment shows that the analysis field generated by NAQPMS-PDAF v2.0 is in good agreement with observations, especially in reducing the underestimation of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ and the overestimation of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC. In particular, the Pearson correlation coefficient (CORR) values for NO${}_{\mathrm{3}}^{-}$, OC, and EC are above 0.96, and the R² values are above 0.93. NAQPMS-PDAF v2.0 also demonstrates superior spatiotemporal interpretation, with most DA sites showing improvements of over 50 %–200 % in CORR and over 50 %–90 % in RMSE for the five chemical components. Compared to the poor performance in the global reanalysis dataset (CORR: 0.42–0.55, RMSE: 4.51–12.27 µg m⁻³) and NAQPMS-PDAF v1.0 (CORR: 0.35–0.98, RMSE: 2.46–15.50 µg m⁻³), NAQPMS-PDAF v2.0 has the highest CORR of 0.86–0.99 and the lowest RMSE of 0.14–3.18 µg m⁻³. The uncertainties in ensemble DA are also examined, further highlighting the potential of NAQPMS-PDAF v2.0 for advancing aerosol chemical component studies.

1 Introduction

PM_2.5 is a complex mixture of various chemical fractions, mainly including sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), nitrate (NO${}_{\mathrm{3}}^{-}$), ammonium (NH${}_{\mathrm{4}}^{+}$), organic carbon (OC), and elemental carbon (EC). These chemical components exert diverse influences on the atmospheric environment (Khanna et al., 2018), human health (Bell et al., 2007; Schlesinger, 2007; J. Li et al., 2022; Alves et al., 2023), and climate change (Schult et al., 1997; Park et al., 2014; Wilcox et al., 2016). However, current detection technologies, such as direct observation by sampling and chemical analysis (Zhang et al., 2015; Ming et al., 2017), ground-based remote sensing inversion (Nishizawa et al., 2008, 2011, 2017), and observation-based machine learning (Lin et al., 2022; Lee et al., 2023; Li et al., 2025), are insufficient in interpreting spatiotemporally continuous information of PM_2.5 chemical components due to the limited number of observation sites or platforms. Although atmospheric chemistry transport models (CTMs) (Wang et al., 2014, 2015; Jia et al., 2017; Yang et al., 2019; Li et al., 2020; Lv et al., 2020) are widely used to characterize the spatiotemporal distribution of multiple chemical species, they are constrained by uncertainties in initial-boundary conditions, physiochemical mechanisms, emission inventories, and meteorological fields (Sax and Isakov, 2003; Mallet and Sportisse, 2006; Rodriguez et al., 2007; Chang et al., 2015; Miao et al., 2020; Xie et al., 2022), resulting in notable discrepancies between the model simulations and accurate observations.

Data assimilation (DA) offers a solution to integrate the multi-source observations, CTMs, and their uncertainties effectively to enhance the simulation and forecasting capabilities of CTMs. Variational methods (3D-Var/4D-Var) (Talagrand and Courtier, 1987), an ensemble Kalman filter (EnKF) (Evensen, 1994, 2003), EnKF variants (EnKFs) (Bishop et al., 2001; Tippett et al., 2003; Hunt et al., 2007; Nerger et al., 2012), and hybrid EnKF–Var methods (Hamill and Snyder, 2000; Schwartz et al., 2014) are most widely applied in DA. However, variational methods have a flow-independent background error covariance (BEC) with the assumption of isotropic, static, and uniform characteristics, and they need to develop a tangent linear adjoint model, which is difficult to practice for complex models. Although EnKFs and hybrid EnKF–Var methods have a flow-dependent BEC, they are sensitive to inadequate ensemble sampling and have high computational costs. Importantly, these methods cannot address model nonlinearity and non-Gaussian error distribution, yielding suboptimal results for DA in highly nonlinear CTMs.

Currently, nonlinear filters, such as the particle filter (PF) (Gordon et al., 1993) and nonlinear ensemble transform filter (NETF) (Tödter and Ahrens, 2015), have been proposed to approximate the complete posterior probability distribution of model states and provide a better representation of non-Gaussian information based on Monte Carlo random sampling and Bayesian theory. However, PF is unstable and susceptible to filter degeneration compared to EnKFs. In a recent study, Nerger (2022) proposed the hybrid Kalman–nonlinear ensemble transform filter (KNETF) to achieve excellent DA performance in the Lorenz-63 and Lorenz-96 models with a smaller ensemble size, which combines the stability of EnKFs and the nonlinearity of NETF (Nerger, 2022). However, to the authors' knowledge, this algorithm has not been applied to the chemical DA of CTMs.

Studies on chemical DA involve the assimilation of aerosol optical properties, such as aerosol optical depth (AOD) and extinction coefficient (EXT), and particulate matter (PM), such as the mass concentrations of PM_2.5 and PM₁₀. The common AOD observations for DA include the OMI AOD (Ali et al., 2013), MODIS AOD (Zhang et al., 2008; Huneeus et al., 2012, 2013; Rubin and Collins, 2014; Lynch et al., 2016; Werner et al., 2019; Kumar et al., 2020), AERONET AOD (Schutgens et al., 2010; Li et al., 2016), sun–sky photometer multiband AOD (Chang et al., 2021), GOCI AOD (Saide et al., 2014; Luo et al., 2020; Kim et al., 2021), and Fengyun/Himawari-8 AOD (Bao et al., 2019; Jin et al., 2019; Xia et al., 2019, 2020). These studies indicated that AOD observations can enhance the accuracy of aerosol simulation and forecast. Compared to AOD, EXT DA effectively improves the interpretation of aerosol vertical distribution (Zhang et al., 2014; Cheng et al., 2019; Wang et al., 2022). Additionally, the simultaneous DA of aerosol optical properties and PM is widely applied in aerosol studies (Tang et al., 2015; Chai et al., 2017). According to our literature review (Yang et al., 2023), there is currently no DA study on aerosol chemical components due to the limited DA influence of PM and AOD on chemical compositions (Chang et al., 2021) and the limited chemical observations with an extensive spatial range. Moreover, the aerosol chemical components exhibit nonlinearity and a non-Gaussian distribution (Ha, 2022), while current mainstream algorithms, such as variational methods or EnKFs, are suboptimal for chemical component DA.

In our previous work, we developed an aerosol vertical DA system (NAQPMS-PDAF v1.0) based on EnKFs to improve the simulation of the extinction coefficient vertical profile (Wang et al., 2022). In this study, we present a novel hybrid nonlinear DA system (NAQPMS-PDAF v2.0) to interpret various PM_2.5 chemical components through online integration of the Parallel Data Assimilation Framework (PDAF; version 2.1, released on 21 February 2023), Observation Module Infrastructure (OMI), and Nested Air Quality Prediction Model System (NAQPMS). We collected 1-month hourly surface observations of five PM_2.5 chemical components (NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, and EC) over northern China and the surrounding areas. We utilized the hybrid localized Kalman–nonlinear ensemble transform filter (LKNETF) to generate a high-resolution and high-accuracy reanalysis dataset of PM_2.5 chemical components for the first time. Notably, the ensemble members in NAQPMS-PDAF v2.0 are generated by perturbing emission species based on their uncertainties and non-Gaussian distribution assumption. Section 2 briefly introduces NAQPMS and PDAF v2.1 with OMI and details the development of NAQPMS-PDAF v2.0, including the system structure, configuration, ensemble generation, and LKNETF algorithm. The data used in this study and experimental settings are also described in Sect. 2. Section 3 presents the DA results, including an evaluation of dependencies, performance, and external comparisons, as well as a discussion of the ensemble DA uncertainty. Section 4 summarizes the conclusions and outlook.

2 Method and data

2.1 NAQPMS

The Nested Air Quality Prediction Modeling System (NAQPMS), developed by the Institute of Atmospheric Physics (IAP), Chinese Academy of Sciences (CAS), is used to provide background fields for key aerosol chemical components in this study. NAQPMS is a multi-scale gridded 3D Eulerian chemical transport model based on continuity equations. The nested grids in the horizontal direction enable data exchange between different domains. Applying terrain-following coordinates in the vertical direction mitigates numerical calculation errors to enhance model accuracy. NAQPMS comprises an input section, a numerical computation section, and an output section. The input section incorporates static terrain data, emission inventories, meteorological fields, and initial-boundary conditions. The numerical computation section performs multiple physicochemical process calculations, including the advection process, eddy diffusion, dry deposition, wet scavenging, gas-phase chemistry, aqueous chemistry, aerosol physicochemical processes (including heterogeneous reactions at the aerosol surface), and other processes. The schemes and features of the physicochemical processes are summarized in Table S1 in the Supplement. The output section is responsible for model post-processing, data diagnostics, and source identification.

NAQPMS is capable of characterizing the 3D spatiotemporal distribution of various atmospheric compositions at global and regional scales and has been widely used in atmospheric pollution and chemistry research, such as O₃ pollution (Wang et al., 2001), haze episodes (Wang et al., 2014; Du et al., 2021), regional transport (Wang et al., 2017, 2019), source identification (Y. Li et al., 2022), air quality simulation at a global scale (Ye et al., 2021) and an urban-street scale (Wang et al., 2023), and acid deposition (Ge et al., 2014).

2.2 PDAF v2.1 with OMI

The Parallel Data Assimilation Framework (PDAF; https://pdaf.awi.de/trac/wiki, last access: 8 March 2024) is an open-source and high-expandability software developed by the Alfred Wegener Institute (AWI) in Germany to integrate observations, numerical models, and assimilation systems for DA tasks and is widely applied in numerical models of meteorology, ocean, land surface, and atmospheric chemistry (Kurtz et al., 2016; Nerger et al., 2020; Mingari et al., 2022; Strebel et al., 2022; Wang et al., 2022; Yu et al., 2022). The initial version of PDAF (PDAF v1.0) was released in 2004. It has undergone continuous improvements and updates, with major updates including the introduction of the ensemble transform Kalman filter (ETKF) and its localized variant (LETKF) in version 1.6; the implementation of PDAF-OMI (Observation Module Infrastructure) in version 1.16; the integration of 3D-Var methods in version 2.0; and the incorporation of the hybrid KNETF and its localized variant (LKNETF) for the first time in version 2.1, which was released in 2023 to handle complex DA situations, such as nonlinearity of the system and non-Gaussian error distribution of the model state. Notably, the version of PDAF coupled in NAQPMS-PDAF v1.0 is PDAF v1.15 (released in 2019), implying that NAQPMS-PDAF v1.0 has more limited applicability and functionality. In this work, PDAF v2.1 is coupled in NAQPMS-PDAF v2.0.

PDAF has offline and online modes. For the offline mode, PDAF and the model perform separately without coupling, obviating the need to modify the model code. For the online mode, PDAF is coupled with the model, and model calculation and data assimilation are performed continuously. Compared to the offline mode, the online coupling has several advantages. Firstly, the initialization of PDAF and the model is integrated, necessitating a single execution rather than two separate executions. Secondly, the model integration result can be directly passed to PDAF for data assimilation. Additionally, the assimilation result of PDAF can be directly passed to the model for the next model integration. The online mode eliminates the need for intermediate steps and improves efficiency. Thirdly, the online mode is controlled by a main program, which allows for efficient use of several processors in the high-performance computing cluster. Conversely, in the offline mode, PDAF and the model are managed by distinct programs, often with fewer processors available for each program. Therefore, the online-mode PDAF is used in this study.

PDAF-OMI, an extension of PDAF, provides I/O interfaces for multi-type observations, simplifying user observation handling by offering generic PDAF-OMI core routines and independent user-supplied routines for each observational type. The user-supplied routines, namely init_dim_obs, init_dim_obs_l, obs_op, and localize_covar, are responsible for reading and writing multi-type observations, applying corresponding observation operators, and performing covariance localization, respectively. The modules for all observation types are integrated into callback_obs_pdafomi, allowing free combinations between different observation types without interference and facilitating collaborative DA for various aerosol chemical components. PDAF-OMI was not applied in NAQPMS-PDAF v1.0. Consequently, NAQPMS-PDAF v1.0 cannot switch between different observational type combinations, and users need to define complete routines for each observation type for the DA process, resulting in more tedious code writing and higher computational costs in NAQPMS-PDAF v1.0.

2.3 NAQPMS-PDAF v2.0

2.3.1 Structure of NAQPMS-PDAF v2.0

Figure 1 illustrates the structure (left portion) and main workflow (right portion) of NAQPMS-PDAF v2.0. As described in the left portion of Fig. 1, the observation part involves the integration of multi-type observations (the purple cuboids) and the utilization of PDAF-OMI. PDAF-OMI enables the simultaneous access and scheduling of multi-type and multi-source observations by employing observational indices, thereby facilitating flexible combinations of observations. The ensemble initial fields (the dark blue cuboids) are crucial inputs for the numerical simulation of NAQPMS. The ensemble forecast/background fields (the dark yellow cuboids) are generated by perturbing emission species based on hypothesized distributions (see Sect. 2.3.3) and performing physiochemical calculations in NAQPMS (the green rectangles). Then, chemical DA is performed by a novel hybrid localized nonlinear DA algorithm (LKNETF; see Sect. 2.3.4) with an adaptive hybrid weight and an adaptive forgetting factor to generate analysis fields (the orange cuboids) for the next realization.

Figure 1 The structure of NAQPMS-PDAF v2.0. (left) The purple cuboids represent the multi-type observations, the dark blue cuboids represent the initial fields, the dark yellow cuboids represent the forecast or background fields, and the orange cuboids represent the analysis fields. Ens.1st represents the first ensemble member, and Ens.Nth represents the Nth ensemble member. (right) The main workflow in NAQPMS-PDAF v2.0, where the blue rectangles represent the modules in NAQPMS, and the yellow rectangles represent the modules in PDAF).

NAQPMS-PDAF v2.0 implements an online coupling between NAQPMS and PDAF v2.1 with OMI, utilizing a level-2 parallel computational framework. The level-2 parallel implementation has been described in our previous work (Wang et al., 2022). The online coupling ensures the continuous operation of model forecasts and assimilation analysis at each time step, achieved by directly integrating PDAF routines into the prototype code of NAQPMS (the right portion of Fig. 1; the blue represents NAQPMS main routines, while the yellow represents PDAF main routines). The level-2 parallel computational framework, which utilizes the message passing interface (MPI) standard, facilitates concurrent processing and data exchange among multiple ensemble members and parallel computation among model state matrixes within each ensemble member, enhancing the efficiency of ensemble analysis and numerical model computations. For instance, the operation of 20 ensemble members necessitates the execution of 20 model tasks, each of which performs integral calculations on a large model grid. A total of 20 model tasks can be executed simultaneously at 20 computational nodes with sufficient computational resources. Each model task can then perform parallel computation with multiple processors by splitting the large model grid into multiple sub-grids. As illustrated in the right portion of Fig. 1, the workflow of NAQPMS-PDAF v2.0 is outlined as follows:

• Step 1. The init_system module initializes NAQPMS by defining all model state variables; allocating numerical matrixes; and configuring parameters, the I/O of meteorological fields, and emission input.

• Step 2. The init_parallel module initializes MPI (MPI_ COMM_WORLD) and the model communicator (MPI_COMM_MODEL), their number of processes, and the rank of a process, followed by init_parallel_pdaf, which initializes MPI communicators for the model tasks, the filter tasks, and the coupling between model and filter tasks.

• Step 3. The initialize module initializes the parameters of the target field, including spatiotemporal dimensions (longitude, latitude, and time steps) and variable dimensions.

• Step 4. The init_pdaf module initializes PDAF variables, such as the local state dimension, global state dimension, and settings for analysis steps.

• Step 5. In this step, the time loop of the forecast and the analysis are performed. The convert_field module is employed to match the matrix storage rule of the target field between NAQPMS and PDAF to ensure compatibility. The field2var module collects the analysis field/initial field and establishes a relationship between the initial field/analysis field and sub-variables in NAQPMS. Subsequently, the analysis field values are allocated to the corresponding NAQPMS sub-variables, and then the NAQPMS_processes module performs the forecast. Afterwards, the var2field module, the inverse of the field2var module, assigns the NAQPMS sub-variables to the forecast field/background field. Finally, the assimilate_pdaf module assimilates the target field with observations to generate an analysis field for the next iteration.

• Step 6. The post-processing module is responsible for finalizing NAQPMS-PDAF, data analysis, and DA evaluation.

2.3.2 Configurations

The meteorological field for NAQPMS is provided by the Weather Research and Forecasting model version 4.0 (WRFv4.0; https://www.mmm.ucar.edu/models/wrf, last access: 26 March 2023). The initial-boundary conditions for WRF are obtained from the NCEP Global Data Assimilation System (GDAS) Final analysis (https://rda.ucar.edu/datasets/ds083.3/, last access: 26 March 2023), with a horizontal resolution of 0.25° × 0.25° and a temporal resolution of 6 h, produced by GDAS. The land use data for WRF were updated by USGS's MCD12Q1 v006 in 2019 (https://lpdaac.usgs.gov/products/mcd12q1v006/, last access: 14 January 2022) with 20 categories. Three nested model domains are conducted with horizontal resolutions of 45 km in the East Asia region (domain 1), 15 km in most areas of China except for the western area (domain 2), and 5 km in the northern China region (domain 3, target research region). WRF and NAQPMS have 40 vertical layers with 27 layers within 2 km. The parameterization schemes for physical processes in WRF are shown in Table S2. The boundary condition input for NAQPMS is provided by the Model for OZone And Related chemical Tracers version 2.4 (MOZART v2.4) global chemistry transport model (Horowitz et al., 2003). The anthropogenic emissions for NAQPMS are from Tsinghua University's 2016 Multi-resolution Emission Inventory for China (MEIC; http://www.meicmodel.org/, last access: 11 March 2023) with a spatial resolution of 0.25° × 0.25°, including residential sources, transportation sources, agricultural sources, industrial sources, and power plant sources. The computational platform is the high-performance supercomputer subsystem cluster with 320 computation nodes, a total of 12 800 processors, and about 153 TB memory on the Big Data Cloud Service Infrastructure Platform (BDCSIP), which meets the demand of NAQPMS-PDAF v2.0 for high-performance parallel computing.

The model state variables include NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, EC, Na⁺, brown carbon, soil PM_2.5, soil PM₁₀, sea salt, fine dust, coarse dust, SO₂, NO₂, and RH. As shown in Fig. 2, the model state has a 4D structure with a longitudinal dimension (ix, 300 grids), latitudinal dimension (iy, 249 grids), variable dimension (ivar, 15), and vertical dimension (iz, 40 layers), in that order. The 4D model state with 15 variables is converted to a 2D state matrix in PDAF; the number of grids in the horizontal-axis direction is ix, and the number of grids in the vertical-axis direction is $iy\times ivar\times iz$. Notably, the 2D state matrix coordinate index contains 3D information for each variable to implement the horizontal and vertical domain localization separately because the horizontal and vertical resolutions are not uniform. This structure has two advantages. First, the parallel cutting of the horizontal axis enables the local domain to retain the full dimensional information ($ix\mathit{\_}p\times iy\times ivar\times iz$, where ix_p is the longitudinal dimension of the local domain). Second, the localization in the local domain permits the analysis to be executed only within a small domain (ix_p×iy) when the length of the horizontal localization radius (R_s) is smaller than iy, effectively reducing the influence of spurious correlations between different state variables. In this study, we set the horizontal and vertical domain localization radius to 200 km (40 grids) and one layer. Additionally, we further implement the observation localization to consider the influence of distance between the analysis grid and observational grid (see Sect. 2.3.4). To minimize computational complexity, the observation errors were assumed to be spatially isotropic, with 0.40, 1.00, 0.50, 3.00, and 0.50 µg m⁻³ for NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, and EC, respectively.

Figure 2 The structure of state variables in NAQPMS-PDAF v2.0.

2.3.3 Generation of ensemble members

In ensemble DA, ensemble members interpret the uncertainty in the model or system, characterized by BEC, which significantly impacts the DA performance (Dai et al., 2014). For CTMs, emission input directly influences the chemical calculation and substantially contributes to the uncertainty. Perturbing emission input can effectively represent the uncertainty in aerosol emissions and enhance the consistency of ensemble error spread, thereby improving aerosol DA (Huang et al., 2023). CTMs are nonlinear, and model state errors are non-Gaussian distributions. To obtain non-Gaussian error distributions, we followed the Kong et al. (2021) method to assume that the emission errors are spatially correlated by an isotropic correlation model with a decorrelation length of 150 km and to generate perturbation coefficient matrixes with the same Gaussian distribution as the emission species, which are subsequently transformed into non-Gaussian distribution matrixes through non-Gaussian process generation v1.2 (Cheynet, 2023).

The target PM_2.5 chemical components are NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, and EC. The perturbed emission species that can directly or indirectly affect the component concentration calculations include SO₂, NO_x, VOC_s, NH₃, CO, PM₁₀, PM_2.5, EC, and OC, with the corresponding uncertainties (δ) listed in Table 1. As shown in Eq. (1), the original emission input matrix (E_p) is multiplied by the corresponding perturbation coefficient matrix (θ_i) to generate the perturbed emission input matrix (E_i) for each emission species. The calculation of the perturbation coefficient matrix (θ_i) is followed by Eqs. (2)–(3). Firstly, N 2D pseudorandom perturbation fields (P_i) are created using Evensen's method (Evensen, 1994). The uncertainties (δ) in the emission species are incorporated into the 2D pseudorandom perturbation fields (P_i) to obtain the final perturbation coefficient matrixes (θ_i). Finally, the Gaussian distribution perturbation coefficient matrixes (θ_i) were transformed into non-Gaussian distribution coefficient matrixes with a given target skewness (set to 1) and kurtosis (set to 6) by non-Gaussian process generation v1.2, which employs the moment-based Hermite transformation model and a cubic transformation.

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathbf{E}}_i={\mathbf{E}}_{\mathrm{p}}\times {\mathit{\theta }}_i,i=\mathrm{1},\mathrm{2},\mathrm{\dots },N\text{(2)}& {\displaystyle }\begin{array}{rl}\ln {\mathit{\theta }}_{o_i}& =({\displaystyle \frac{({\mathbf{P}}_i-\frac{\mathrm{1}}{N}\times {\sum }_{i=\mathrm{1}}^N{\mathbf{P}}_i)}{\sqrt{\frac{\mathrm{1}}{N}\times {\sum }_{i=\mathrm{1}}^N({\mathbf{P}}_i-\frac{\mathrm{1}}{N}\times {\sum }_{i=\mathrm{1}}^N{\mathbf{P}}_i{)}^{\mathrm{2}}}}}\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\times \ln (\mathrm{1}+{\mathit{\delta }}^{\mathrm{2}}))\times \sqrt{\ln (\mathrm{1}+{\mathit{\delta }}^{\mathrm{2}})}\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathit{\theta }}_i& ={\displaystyle \frac{\left({\mathit{\theta }}_{o_i}-\frac{\mathrm{1}}{N}\times {\sum }_{i=\mathrm{1}}^N{\mathit{\theta }}_{o_i}\right)}{\sqrt{\frac{\mathrm{1}}{N}\times {\sum }_{i=\mathrm{1}}^N{\left({\mathit{\theta }}_{o_i}-\frac{\mathrm{1}}{N}\times {\sum }_{i=\mathrm{1}}^N{\mathit{\theta }}_{o_i}\right)}^{\mathrm{2}}}}}\\ & \times \left({\displaystyle \frac{\mathrm{1}}{N}}\times {\sum }_{i=\mathrm{1}}^N{\mathit{\theta }}_{o_i}\right)\times \mathit{\delta }+{\displaystyle \frac{\mathrm{1}}{N}}\times {\sum }_{i=\mathrm{1}}^N{\mathit{\theta }}_{o_i}\end{array}\end{array}$$

Notably, all matrix operations involved are Schur products. Here, E_i denotes the ith ensemble perturbed emission input matrix, E_p indicates the original unperturbed emission input matrix, and θ_i represents the ith ensemble perturbation coefficient matrix. ${\mathit{\theta }}_{o_i}$ is the ith ensemble original perturbation coefficient matrix, which is obtained by mathematical transformation of the ith ensemble pseudorandom perturbation matrix P_i, including standardization and scaling by uncertainty (δ), as well as a logarithm.

Table 1 The uncertainties in emission species in NAQPMS-PDAF v2.0.

2.3.4 Hybrid nonlinear DA algorithm with adaptive forgetting factor

To thoroughly integrate the stability of EnKFs with the nonlinearity of nonlinear filters and be ideal for the nonlinear and non-Gaussian distribution situations, the hybrid LKNETF is used in this study. This section reviews the algorithms of LETKF, LNETF, and their combination (LKNETF).

ETKF, a deterministic filter in EnKFs, efficiently obtains analysis samples using a transformation matrix and the square root of the forecast error covariance (Bishop et al., 2001). In contrast to stochastic filters in EnKFs, ETKF prevents underestimation of the analysis error covariance resulting from the random observation perturbations. And it is particularly applicable in situations with small ensemble sizes (Lawson and Hansen, 2004). The realization of ETKF can be divided into the forecast and analysis steps.

In the forecast step, the forecast state vector (${\mathit{x}}_t^{\mathrm{f}}$) at t is generated by numerical model (M) integration of the analysis state vector (${\mathit{x}}_{t-\mathrm{1}}^{\mathrm{a}}$) at t−1. The forecast error covariance matrix (${\mathbf{P}}_t^{\mathrm{f}}$) can be calculated by the perturbation of the forecast ensemble (${\mathit{X}}_t^{{\mathrm{f}}^{\prime }}$):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{x}}_t^{\mathrm{f}}=\mathbf{M}\left({\mathit{x}}_{t-\mathrm{1}}^{\mathrm{a}}\right),{\mathit{X}}_t^{\mathrm{f}}=[{\mathit{x}}_{{\mathrm{1}}_t}^{\mathrm{f}},{\mathit{x}}_{{\mathrm{2}}_t}^{\mathrm{f}},\mathrm{\dots },{\mathit{x}}_{K_t}^{\mathrm{f}}],\text{(5)}& {\displaystyle }{\mathbf{P}}_t^{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{(K-\mathrm{1})}}{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}{{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}}^{\mathrm{T}},\end{array}$$

where ${\mathit{X}}_t^{\mathrm{f}}$ is the forecast ensemble at t, and K is the number of ensemble members. ${\mathit{X}}_t^{{\mathrm{f}}^{\prime }}$ is the perturbation of the forecast ensemble at t, calculated by ${\mathit{X}}_t^{\mathrm{f}}$, and the forecast ensemble mean $\stackrel{\mathrm{‾}}{{\mathit{X}}_t^{\mathrm{f}}}$ at t.

In the analysis step, the forecast error covariance matrix (${\mathbf{P}}_t^{\mathrm{f}}$) at t is transformed to the analysis error covariance matrix (${\mathbf{P}}_t^{\mathrm{a}}$) at t by a transform matrix (T):

$$\begin{array}{}\text{(6)}& {\mathbf{P}}_t^{\mathrm{a}}={\mathit{X}}_t^{{\mathrm{f}}^{\prime }}\mathbf{T}{{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}}^{\mathrm{T}}.\end{array}$$

The transform matrix (T) is defined as follows and can be decomposed to a left singular vector matrix (U), a singular value matrix (S), and a right singular vector matrix (V) through the singular value decomposition:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathbf{T}}^{-\mathrm{1}}={\mathit{\rho }}_{\mathrm{adaptive}}\left(K-\mathrm{1}\right)\mathbf{I}+\left(\mathbf{H}{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}{)}^{\mathrm{T}}\right(\mathbf{L}\cdot {\mathbf{R}}^{-\mathrm{1}})\mathbf{H}{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}=\mathbf{USV},\text{(8)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{adaptive}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{ens}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{resid}}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{obs}}^{\mathrm{2}}}},\end{array}$$

where ρ_adaptive is an adaptive forgetting factor used for the inflation of error covariance estimation (the initial ρ_adaptive is set to 0.9 in this study). ${\mathit{\sigma }}_{\mathrm{ens}}^{\mathrm{2}}$ is the mean ensemble variance, ${\mathit{\sigma }}_{\mathrm{resid}}^{\mathrm{2}}$ is the mean of the observation-minus-forecast residual, and ${\mathit{\sigma }}_{\mathrm{obs}}^{\mathrm{2}}$ is the mean observation variance. I is the identity matrix. H is the observation operator. L is the localization matrix, a weight matrix calculated by the fifth-order polynomial (Nerger, 2015), implemented in LETKF for observation localization analysis to avoid observational spurious correlation and filter divergence effectively (Hunt et al., 2007). R is the observation error covariance matrix.

The analysis state vector (${\mathit{x}}_t^{\mathrm{a}}$) at t is calculated by the forecast state vector (${\mathit{x}}_t^{\mathrm{f}}$) at t, the perturbation of the forecast ensemble (${\mathit{X}}_t^{{\mathrm{f}}^{\prime }}$) at t, and a weight vector (w):

$$\begin{array}{}\text{(9)}& {\mathit{x}}_t^{\mathrm{a}}={\mathit{x}}_t^{\mathrm{f}}+{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}\mathit{w}.\end{array}$$

The weight vector (w) is given by the following equation:

$$\begin{array}{}\text{(10)}& \mathit{w}=\mathbf{T}\left(\mathbf{H}{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}{)}^{\mathrm{T}}\right(\mathbf{L}\cdot {\mathbf{R}}^{-\mathrm{1}}\left)\right(\mathit{y}-\mathbf{H}{\mathit{x}}_t^{\mathrm{f}}),\end{array}$$

where y is observations.

The analysis ensemble (${\mathit{X}}_t^{\mathrm{a}}$) at t can be obtained by the analysis ensemble mean ($\stackrel{\mathrm{‾}}{{\mathit{X}}_t^{\mathrm{a}}}$) at t, the perturbation of the forecast ensemble (${\mathit{X}}_t^{{\mathrm{f}}^{\prime }}$) at t, and a transform matrix (C) represented by the symmetric square root of T:

$$\begin{array}{}\text{(11)}& {\mathit{X}}_t^{\mathrm{a}}=\stackrel{\mathrm{‾}}{{\mathit{X}}_t^{\mathrm{a}}}+\sqrt{K-\mathrm{1}}{\mathit{X}}_t^{{\mathrm{f}}^{\prime }}\mathbf{C}.\end{array}$$

The transform matrix (C) is calculated as follows:

$$\begin{array}{}\text{(12)}& \mathbf{C}={\mathbf{US}}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{U}}^{\mathrm{T}}.\end{array}$$

NETF is a second-order exact ensemble square root filter effectively applied to the nonlinear and non-Gaussian DA (Tödter and Ahrens, 2015). Like PF, NETF indirectly updates the model state by using observations to affect the weights of the prior ensemble. However, PF and NETF differ in the sampling method. PF utilizes the Monte Carlo and Bayesian approaches to calculate particle weights based on observations, which are then used to generate the analysis ensemble by weighting the resampling forecast ensemble. In high-dimensional systems, as the DA progresses, the weight differences among particles increase, with most particles having weights close to 0, leading to filter degeneration. In contrast, NETF generates the analysis ensemble through a deterministic matrix square root transformation of the forecast ensemble, where the mean and covariance matrix of the analysis ensemble match the weighted values in PF (as shown in Eqs. 13–14). Due to the similarity between NETF and ETKF, the localization can be implemented in NETF (LNETF) (Tödter et al., 2016).

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{K}}{\sum }_{i=\mathrm{1}}^K{\mathit{x}}_i^{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{K}}{\sum }_{i=\mathrm{1}}^K{\mathit{w}}_i{\mathit{x}}_i^{\mathrm{f}}\end{array}$$

Here, ${\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{a}}$ is the analysis state vector mean, K is the number of ensemble members, ${\mathit{x}}_i^{\mathrm{a}}$ is the ith analysis state vector, w_i is the ith particle weight in PF (which is calculated by the Bayesian method ${\mathit{w}}_i=\frac{p\left(\mathit{y}|{\mathit{x}}_i^{\mathrm{f}}\right)}{p\left(\mathbf{y}\right)}$), y is the observations, and ${\mathit{x}}_i^{\mathrm{f}}$ is the ith forecast state vector.

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathbf{P}}^{\mathrm{a}}& ={\displaystyle \frac{\mathrm{1}}{K-\mathrm{1}}}{\sum }_{i=\mathrm{1}}^K({\mathit{x}}_i^{\mathrm{a}}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{a}})({\mathit{x}}_i^{\mathrm{a}}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{a}}{)}^{\mathrm{T}}\\ & ={\sum }_{\mathrm{i}=\mathrm{K}}^K{\mathit{w}}_i({\mathit{x}}_i^{\mathrm{f}}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{f}})({\mathit{x}}_i^{\mathrm{f}}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{f}}{)}^{\mathrm{T}}\end{array}\end{array}$$

Here, P^a is the error covariance matrix of the analysis ensemble, calculated by the perturbation of the analysis ensemble.

In NETF, A acts as a transform matrix like the transform matrix (T) in ETKF, which can be obtained from the weight vector (w).

$$\begin{array}{}\text{(15)}& {\displaystyle }{\mathbf{P}}^{\mathrm{a}}={\displaystyle \frac{K}{(K-\mathrm{1})}}{\mathbf{X}}^{{\mathrm{f}}^{\prime }}{{\mathbf{AX}}^{{\mathrm{f}}^{\prime }}}^{\mathrm{T}}\text{(16)}& {\displaystyle }{\mathbf{A}}^{\frac{\mathrm{1}}{\mathrm{2}}}={\left(\mathbf{W}-\mathit{w}{\mathit{w}}^{\mathrm{T}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}={\mathbf{VD}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}^{\mathrm{T}}\end{array}$$

Here, the matrix W≡diag(w) is defined as a diagonal matrix created from the weight vector (w). A can be decomposed (A=VDV^T) by a singular value decomposition as it is a real, symmetric, positive semidefinite matrix. V is the orthogonal matrix, and D is a diagonal matrix.

Then, the perturbation of the analysis ensemble (${\mathbf{X}}^{{\mathrm{a}}^{\prime }}$) and the analysis ensemble (X^a) can be obtained by applying the square root of A as a transform matrix:

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathbf{X}}^{{\mathrm{a}}^{\prime }}=\sqrt{K}{\mathbf{X}}^{{\mathrm{f}}^{\prime }}{\mathbf{A}}^{\frac{\mathrm{1}}{\mathrm{2}}},\text{(18)}& {\displaystyle }{\mathbf{X}}^{\mathrm{a}}={\stackrel{\mathrm{‾}}{\mathbf{X}}}^{\mathrm{f}}+{\mathbf{X}}^{{\mathrm{f}}^{\prime }}\left(\stackrel{\mathrm{‾}}{\mathbf{W}}+\sqrt{K}{\mathbf{A}}^{\frac{\mathrm{1}}{\mathrm{2}}}\right).\end{array}$$

LKNETF combines LETKF and LNETF through a hybrid weight γ to perform better in systems with different nonlinearity degrees and to be implemented in situations with smaller ensemble sizes (Nerger, 2022). When γ approaches 1, the analysis increment (ΔX_LETKF) computed by LETKF becomes more significant and appropriate for linear systems with Gaussian distributions. Conversely, when γ approaches 0, the analysis increment (ΔX_LNETF) computed by LNETF becomes more significant and appropriate for nonlinear systems with non-Gaussian distributions. The one-step update scheme is used in this study.

$$\begin{array}{}\text{(19)}& {\mathbf{X}}_{\mathrm{HSync}}^{\mathrm{a}}={\stackrel{\mathrm{‾}}{\mathbf{X}}}^{\mathrm{f}}+\left(\mathrm{1}-\mathit{\gamma }\right)\mathrm{\Delta }{\mathbf{X}}_{\mathrm{LNETF}}+\mathit{\gamma }\mathrm{\Delta }{\mathbf{X}}_{\mathrm{LETKF}}\end{array}$$

2.4 Data

2.4.1 Observation

A total of 1 month's (February 2022) hourly mass concentration observations of five PM_2.5 chemical components (NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, and EC) from 33 ground-based sites in northern China and the surrounding areas were collected for this work (Fig. 3). Out of the 33 sites, 24 (DA sites) were utilized for DA and internal validation, and the remaining 9 (VE sites) were used for independent verification to assess the influence of DA sites on neighboring areas. These sites were divided using the K-means clustering algorithm (Lloyd, 1982; Arthur and Vassilvitskii, 2007). The Supplement provides a detailed description (Sect. S1). PM_2.5 hourly observations from the China National Environmental Monitoring Centre (CNEMC; http://www.cnemc.cn/, last access: 1 November 2023) were employed to assess the overall mass concentration of PM_2.5 chemical components in NAQPMS-PDAF v2.0. Due to incomplete spatial overlap between the PM_2.5 sites and the chemical component sites, the PM_2.5 sites were selected based on the closest Euclidean distance between PM_2.5 sites and chemical component sites.

Figure 3 The model domains in the WRF simulation (a) and the location of the observations (b). Domain 3 in (a) is the target area of this study. The 24 red sites in (b) represent the sites for data assimilation, and the 9 green sites in (b) represent the sites for independent spatial validation. The topographic dataset is from the ETOPO1 1 arcmin global relief model, taken from the National Geophysical Data Center (NOAA National Geophysical Data Center, 2009).

2.4.2 Global reanalysis dataset

The global reanalysis datasets of PM_2.5 chemical components in February 2022 were obtained from the Copernicus Atmosphere Monitoring Service ReAnalysis (CAMSRA; 0.75° × 0.75°) (Inness et al., 2019) and the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2; 0.5° × 0.625°) (Randles et al., 2017) to compare with the reanalysis dataset generated by NAQPMS-PDAF v2.0. For the data consistency, the global reanalysis surface grid data located in the observation sites of the PM_2.5 chemical component were extracted through the k-nearest-neighbor search method (Friedman et al., 1977), which can efficiently match grid points and observation sites based on latitude and longitude data and Euclidean distances. Our 3-hourly NAQPMS-PDAF v2.0 output of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ was extracted to compare with the CAMSRA dataset, and hourly NAQPMS-PDAF v2.0 output of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC was extracted to compare with the MERRA-2 M2T1NXAER dataset.

2.5 Experimental setting and evaluation method

In our study, four tests were conducted to evaluate the performance of NAQPMS-PDAF v2.0 with hourly observations of five PM_2.5 chemical components, including (1) the dependence on ensemble size and assimilation frequency, (2) the ability to interpret mass concentration and spatiotemporal characteristics, (3) the quality of output data compared to other reanalysis datasets, and (4) the uncertainty in ensemble assimilation. In practice, the ratio of the ensemble size to the number of processes – 1:50 – in high-performance computers was the optimal parallel scheme to balance computing efficiency and computing resources (Wang et al., 2022).

All the tests were run in NAQPMS-PDAF v2.0 after a spin-up experiment with 24 time steps from 00:00 to 23:00 (LST) on 1 February 2022. (1) For the first test, we assimilated the hourly observations of five PM_2.5 chemical components from all sites with 48 time steps from 00:00 (LST) on 2 February to 23:00 (LST) on 3 February 2022. In the first scenario, we controlled a fixed assimilation frequency of 1 h and changed the ensemble size to 2, 5, 10, 15, 20, 30, 40, and 50. In the second scenario, we controlled a fixed ensemble size of 20 and changed the assimilation frequency to 1, 2, 3, 4, 5, 6, 8, and 12 h. (2) For the second test, we set an ensemble size of 20 and an assimilation frequency of 1 h and assimilated the hourly observations of five PM_2.5 chemical components from DA sites with 648 time steps from 00:00 (LST) on 2 February to 23:00 (LST) on 28 February 2022. We also conducted a free-running (FR) experiment without assimilation in the same period for comparison. (3) For the third test, we followed the settings in the second test but assimilated the observation from all sites to generate a high-quality reanalysis dataset of five PM_2.5 chemical components. (4) The final test was analogous to the first test but with a distinct scenario designed to examine the influence of ensemble perturbation on ensemble assimilation. From Table 2, we fixed species uncertainty (M4 setting) with five distribution types in the first scenario and fixed distribution type (T2 setting) with five SO₂ uncertainties in the second.

Table 2 The experiment settings for emission perturbation.

We used the continuous ranked probability score (CRPS) to evaluate ensemble size dependency, which measures the consistency between the ensemble forecast distribution and corresponding observations (Jolliffe and Stephenson, 2012). The calculation rules are referred to in Hersbach's study (Hersbach, 2000). Moreover, four common statistical indicators, the Pearson correlation coefficient (CORR), root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²), were used to assess the DA system performance in interpreting PM_2.5 chemical components (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, OC, and EC). The CORR measures the correlation between the system outputs and corresponding observations, the RMSE and MAE indicate the overall system accuracy, and the R² reflects the proportion of variability in the observations explained by the assimilation system.

3 Results and discussion

3.1 The dependence on ensemble size and assimilation frequency for five components

Ensemble size is a crucial parameter in ensemble assimilation, determining the model state's uncertainty range. A larger ensemble size more accurately represents the error distribution of state variables but requires considerable computing resources and time, especially for high-dimensional systems. A smaller ensemble size can easily lead to underestimating the error covariance matrix, especially for the fine-resolution model (Kong et al., 2021). Thus, identifying an appropriate ensemble size to balance computational efficiency and accuracy is the primary step in ensemble DA. Our prior study (NAQPMS-PDAF v1.0) only evaluated the correlation between ensemble size and parallel efficiency and concluded that the ratio of ensemble size to high-performance computing processors was 1:50 (Wang et al., 2022), while the impact of ensemble size on the accuracy and computational efficiency was neglected. This study assesses the NAQPMS-PDAF v2.0 dependency on ensemble size through three statistical indicators (CRPS, RMSE, and CORR). Figure 4 shows the mean CRPS, RMSE, and CORR values and the statistical averages of the elapsed time over 48 time steps with the ensemble sizes of 2, 5, 10, 15, 20, 30, 40, and 50.

From Fig. 4a, when the ensemble size is at its minimum level of 2, the mean CRPS values of the five PM_2.5 chemical components are more significant, with NO${}_{\mathrm{3}}^{-}$ exhibiting the most considerable difference between the simulation distribution and observations (more than 4). With each increase in ensemble size, the mean CRPS values of the five chemical components progressively reduce and eventually reach convergence when the ensemble size is 10, implying that a hybrid nonlinear filter can maintain high accuracy and reliability in ensemble assimilation with an ensemble size that is smaller than the traditional minimum of 20 ensemble members, as observed in prior ensemble assimilation studies (Constantinescu et al., 2007; Miyazaki et al., 2012; Schwartz et al., 2014; Rubin et al., 2017; Kong et al., 2021; Tsikerdekis et al., 2021; Wang et al., 2022), including NAQPMS-PDAF v1.0. The mean CRPS value of EC is the lowest among the five chemical components, indicating the highest accuracy and reliability of EC ensemble DA. The performance of other components is similar. Like CRPS values, the values of RMSE and CORR decrease and increase, respectively, as the ensemble size increases, and convergence begins to occur when the ensemble size is 10 (Fig. 4b and c). Compared with other chemical components, the CORR value of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ is significantly lower, less than 0.8, possibly due to its estimated background field error covariance driven by the inadequate ensemble perturbations. Therefore, in the Discussion section, we discuss the uncertainties in ensemble perturbations.

Figure 4 Assessment of ensemble size dependency based on mean continuous ranked probability score (CRPS) (a), root mean square error (RMSE) (b), correlation coefficient (CORR) (c), and time (d).

Figure 4d shows the time required for the four processes of ensemble assimilation under different ensemble sizes, including initialization, model integration, assimilation, and post-processing. The model integration process in NAQPMS-PDAF v2.0 takes the longest, followed by post-processing, initialization, and assimilation. The required time for initialization and post-processing increases with increasing ensemble size, while for model integration and assimilation, except for ensemble size 30, the required time is the same under different ensemble sizes. Generally, the time needed for ensemble sizes of 30–50 is considerably higher than that for smaller ones. Although convergence occurs with an ensemble size of 10, our work illustrates a similar time required between ensemble sizes 10 and 20. Consequently, we selected an ensemble size of 20 to ensure optimal performance of NAQPMS-PDAF v2.0, considering both assimilation efficiency and accuracy.

The assimilation frequency is the interval at which observational data are introduced into the DA system, directly affecting the practical assimilation data volume and computation cost. High-frequency DA with high-quality observations is crucial for improving numerical simulations and forecasts (Liu et al., 2021). Figure 5 demonstrates that the MAE values of the five chemical component analysis fields range from 0.02 to 0.12 µg m⁻³, RMSE values range from 0.23 to 2.61 µg m⁻³, and CORR values range from 0.71 to 0.98 at a 1 h assimilation time interval, which is significantly better than the statistical indicators at lower assimilation frequencies. Even at a 2 h assimilation frequency, the assimilation effect drops sharply compared to the 1 h interval, especially for NO${}_{\mathrm{3}}^{-}$, OC, and EC. The values of MAE and RMSE increase by 2.6–5.82 and 4.72–9.57 µg m⁻³, respectively, and the CORR values decrease by 0.27–0.81. Gradual increasing trends in MAE and RMSE values and a slight decreasing trend in CORR values are observed as the assimilation frequency decreases from the 2 h interval. Therefore, the fast-updating assimilation with a 1 h interval significantly improves the NAQPMS simulation. For the forecasting field (Fig. S2), the low sensitivity of state variables to assimilation frequency suggests that NAQPMS-PDAF v2.0 can appropriately reduce assimilation frequency during the actual forecasting phase, lowering the demand for high-temporal-resolution observations and computational resources.

Figure 5 Assessment of assimilation interval dependency based on mean absolute error (MAE) (a), root mean square error (RMSE) (b), and correlation coefficient (CORR) (c) at the analysis step.

3.2 Evaluation of NAQPMS-PDAF v2.0 performance

3.2.1 Overall validation of DA results

We conducted a control experiment (free-running (FR) field) without any DA and with a DA experiment. This section verifies the forecast (FOR) field and analysis (ANA) field at 24 DA sites and 9 VE sites, respectively. Figure 6 shows the scatter distribution of observations and simulations at DA sites. For the FR field (Fig. 6a1–a5), five chemical components have CORR values ranging from 0.32 to 0.56, and R² values do not exceed 0.3, indicating poor consistency between observations and simulations. In detail, the simulated mass concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC are significantly overestimated, while the simulated concentrations of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ are underestimated. OC has the most significant error, with an RMSE value of 25.84 µg m⁻³ and an MAE value of 19.41 µg m⁻³. Moreover, the error distributions of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ are close to a symmetric distribution with a mean value of 0, while the error distributions of OC and EC are skewed to the left from the mean value of 0 (Fig. 7a1–a5), showing the relatively better simulations in SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ than in OC and EC. Overall, NAQPMS cannot interpret the mass concentrations of the five chemical components with significant errors, mainly due to the uncertainties in chemical mechanisms (Miao et al., 2020).

Figure 6 Scatterplots of the DA site simulations versus the DA site observations with probability density for the free-running (FR) field (a1–a5), forecast (FOR) field (b1–b5), and analysis (ANA) field (c1–c5). The stippled gray lines represent the 2:1, 1:1, and 1:2 lines, and the solid red line represents the fitting regression line.

Figure 7 Probability distributions of bias between DA site observations and DA site simulations for the free-running (FR) field (a1–a5), forecast (FOR) field (b1–b5), and analysis (ANA) field (c1–c5).

After DA, FOR shows a slight improvement with a slight increase in CORR and R² and a decrease in RMSE and MAE, especially for NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ (Fig. 6b1–b5). Although SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC are significantly overestimated with a slight decrease in CORR and R², the RMSE and MAE values decrease. Additionally, the error distributions of the five chemical components are concentrated at 0, and the overestimation of OC and EC has been improved compared to FR (Fig. 7b1–b5). These results indicate that DA reduces the overall FOR errors in NAQPMS due to improved forecasting ability by obtaining optimal initial fields. However, further improvements are necessary to address the NAQPMS uncertainties in emission sources, meteorological input, and imperfect physiochemical mechanisms. For ANA (Fig. 6c1–c5), DA significantly improves the simulations of the five chemical components, making the ANA consistent with the observations. The CORR values are not less than 0.86; the RMSE and MAE values do not exceed 3.23 and 1.49 µg m⁻³, respectively; and the R² values are not less than 0.74. Specifically, the CORR values for NO${}_{\mathrm{3}}^{-}$, OC, and EC are not less than 0.96, and the R² values are not less than 0.93. The error distributions of the five chemical components concentrate to 0, with the mean bias ranging from 0 ± 0.08 to 1.02 ± 3.07 µg m⁻³ (Fig. 7c1–c5). The results of VE sites show similar characteristics to the DA sites (Figs. S3 and S4). Compared to FR, the overall errors of the FOR and ANA for the five chemical components decrease with a significant improvement in ANA, showing that the CORR values of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ increase by 0.15 and 0.45, respectively; the R² values of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ increase by 0.22 and 0.81, respectively; and the RMSE values of OC and EC decrease by 21.77 and 17.79 µg m⁻³, respectively. Overall, the FOR and ANA errors decreased significantly. The ANA of the five chemical components at DA sites is almost entirely consistent with the observations, indicating excellent DA performance.

3.2.2 Assessment of temporal variation in chemical components

The ensemble DA employs a cyclic updating process wherein the forecast and analysis steps are continuously completed at each iteration (Evensen, 2003; Houtekamer and Zhang, 2016). In the forecast step, the ANA at the current time step serves as the optimal initial field to advance the model integration and obtain the FOR at the next step. In the analysis step, the FOR at the next time step provides background field information for the subsequent DA analysis to generate the ANA at the next time step. The FOR and ANA interact with each other in the temporal dimension. Therefore, in this section, we assess the ability of NAQPMS-PDAF v2.0 to interpret the temporal variations in the five chemical components. Figure 8 illustrates the time series of the five chemical components at two representative sites, including a DA site in Tianjin and a VE site in Heze. For the DA site (Fig. 8a), the temporal variations in NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ in FR and FOR exhibit better agreement with the observed temporal variations (OBS) than those of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC. However, NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ mass concentrations are significantly lower than the high-value mass concentrations observed on 25 February. The mass concentration of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ in FR is greatly overestimated during the periods of 9–11, 18–19, and 24–25 February. The mass concentrations of OC and EC in FR are overestimated throughout February with substantial temporal fluctuations. Although the time series of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC in FOR show some improvement, noticeable differences from the OBS are still apparent. After DA, the ANA time series for the five chemical components align well with the OBS, indicating good consistency and accurate representation of temporal characteristics, such as the NH₄NO₃ pollution captured on 25 February. Notably, the mass concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ peaked on 10–11 and 26 February, indicating intensified atmospheric secondary chemical reactions, primarily due to neutralization reactions of acidic pollutants capturing NH₃. The temporal variations in NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ are more similar because atmospheric NO${}_{\mathrm{3}}^{-}$ mainly exists as NH₄NO₃ rather than other metal nitrates, and NH₄NO₃ can form before the complete neutralization of H₂SO₄ (Ge et al., 2017). The improvements at the VE site (Fig. 8b) are like those at the DA site, with the ANA time series of the five chemical components showing closer agreement with the OBS, which suggests that the localization analysis in DA effectively facilitates the propagation of observations within a specific spatial range and mitigates the assimilation anomalies caused by spurious correlations from the distant sites (Hunt et al., 2007).

Figure 8 Hourly variation in five PM_2.5 chemical components in a representative DA site (a) and a representative VE site (b).

NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, and EC are critical chemical components of PM_2.5, and the sum of their mass concentrations can be approximated as the PM_2.5 mass concentration. We further assessed the simulation enhancement of PM_2.5 time series based on ground-level PM_2.5 observations. Six representative sites were selected, including three DA sites (Fig. 9a1–a3) and three VE sites (Fig. 9b1–b3). The FR and FOR in DA and VE sites show significant overestimation and poor consistency with the OBS, mainly due to the overestimation of OC and EC mass concentrations. Conversely, the PM_2.5 time series in ANA closely matches that of the OBS, accurately capturing the actual variation in PM_2.5. In some specific instances, such as on 26 February at 00:00 in Tianjin and Langfang, the peak value of ANA was lower than that of the OBS, which could be attributed to the negligence of other PM_2.5 components (such as mineral dust and sea salt) and the inconsistency in location between ground-level PM_2.5 observational sites and chemical component observational sites. Overall, the DA of chemical component observations significantly enhanced the simulation of PM_2.5 time series in NAQPMS. Compared to the CORR values of FR and FOR, the CORR values of ANA at the six representative sites increased by 13.64 %–89.58 % and 17.19 %–75.00 %, while the RMSE values decreased by 56.03 %–83.13 % and 40.74 %–72.20 % (Table S3).

Figure 9 Hourly variation in PM_2.5 in three representative DA sites (a1–a3) and three representative VE sites (b1–b3).

3.2.3 Assessment of spatial distribution in chemical components

DA can improve the interpretation of model states in the analysis domain by using a limited number of observations. The ability to represent spatial distribution accurately is a crucial function for aerosol DA. Figure 10 displays the spatial distribution of the monthly average mass concentrations for the five chemical components, including OBS, FR, FOR, ANA, and analysis increment (INC). The spatial distributions of bias and statistical indicators for FR, FOR, and ANA are shown in Figs. 11 and 12, respectively.

Figure 10 Spatial concentration distribution of the site observation (OBS) (a1–e1), free-running (FR) field (a2–e2), forecast (FOR) field (a3–e3), analysis (ANA) field (a4–e4), and increment (INC) between ANA and FR (a5–e5) for five PM_2.5 chemical components.

Figure 11 Spatial distribution of DA site bias for five PM_2.5 chemical components from the observation (OBS) for the free-running (FR) field (a1–e1), forecast (FOR) field (a2–e2), and analysis (ANA) field (a3–e3).

Figure 12 Spatial distribution of DA site statistical indicators for five PM_2.5 chemical components. Panels (a1)–(e1) represent the values of RMSE and CORR for the forecast (FOR) field, (a2)–(e2) are the same as (a1)–(e1) but for the analysis (ANA) field, (a3)–(e3) represent the improvement in RMSE and CORR for the forecast (FOR) field, and (a4)–(e4) are the same as (a3)–(e3) but for the analysis (ANA) field. The size represents the value of RMSE in (a1)–(e2) and the improvement percentage compared to non-assimilation in (a3)–(e4).

The spatial characteristics of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ are similar. Compared to the OBS (Fig. 10a1 and c1), the FR (Fig. 10a2 and c2) and FOR (Fig. 10a3 and c3) have failed to capture the high-value mass concentrations in the border area between Hebei Province, Shanxi Province, Henan Province, and Shandong Province, especially in the northern region of Henan Province. The primary reason is the uncertainties in emission inventories in winter heating periods, which result in insufficient emission statistics of gaseous precursors NO_x and NH₃ (Aleksankina et al., 2018). After DA, this situation is significantly improved with the ANA (Fig. 10a4 and c4). The INCs in the Beijing–Tianjin–Hebei region, Shanxi Province, Henan Province, and Shandong Province are positive (Fig. 10a5 and c5), indicating varying degrees of improvement in correcting the underestimation of mass concentrations. Specifically, for NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ at DA sites, the biases between the OBS and ANA are significantly reduced compared to the biases between the OBS and FR (Fig. 11), with the mean absolute bias decreasing by 0.93 and 4.27 µg m⁻³, respectively. Moreover, the overall biases at VE sites also decrease (Fig. S5). As for the spatial statistical indicators of NH${}_{\mathrm{4}}^{+}$ (Fig. 12a1 and a2), the CORR values in FOR and ANA range from 0.39 to 0.79 and 0.70 to 0.97, respectively, and the RMSE values range from 3.16 to 7.65 µg m⁻³ and 1.20 to 3.49 µg m⁻³, respectively. As for the spatial statistical indicators of NO${}_{\mathrm{3}}^{-}$ (Fig. 12c1 and c2), the CORR values in FOR and ANA range from 0.09 to 0.76 and 0.89 to 0.99, respectively, and the RMSE values range from 4.88 to 15.69 µg m⁻³ and 1.34 to 5.39 µg m⁻³, respectively. For the FOR, the improvement in accuracy for NO${}_{\mathrm{3}}^{-}$ is more significant than that for NH${}_{\mathrm{4}}^{+}$, with the CORR values of most DA sites increasing by more than 10 % and the RMSE of most DA sites decreasing by no less than 10 % (Fig. 12a3 and c3). For the ANA, NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ exhibit significant improvements in CORR and RMSE, as most DA sites show over 150 % improvement in CORR and over 50 % improvement in RMSE (Fig. 12a4 and c4). Improvements can also be found in NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ at VE sites (Fig. S6). The spatial consistency of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ indicates that NH₄NO₃ is the primary aerosol chemical component, highlighting the necessity of coordinated control of precursor NO_x and NH₃.

Unlike NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$, compared to the OBS (Fig. 10b1), the mass concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ in the FR and FOR (Fig. 10b2 and b3) are significantly overestimated, especially in Shandong Province. In contrast, the ANA has improved dramatically (Fig. 10b4), with most areas showing negative INCs (Fig. 10b5). The mean absolute biases in DA and VE sites have decreased by 1.80 and 2.68 µg m⁻³, respectively (Figs. 11 and S5). Specifically, after DA, the CORR values of the FOR and ANA range from 0.22 to 0.71 and 0.58 to 0.97, and the RMSE values range from 3.42 to 11.07 µg m⁻³ and 1.20 to 4.30 µg m⁻³, respectively (Fig. 12b1 and b2). The CORR and RMSE values in FOR improved significantly (Fig. 12b3) at DA sites around Beijing, while the CORR values in ANA increased by more than 13 %, with most DA sites showing an increase of over 50 %, and RMSE values decreased by no less than 30 %, with most DA sites showing a decrease of over 70 % (Fig. 12b4). Moreover, half of the VE sites show significant improvement in the CORR and RMSE in the FOR and ANA, mainly due to their proximity to more DA sites (Fig. S6). The OBS and ANA indicate a considerable control in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ pollution during the winter heating period due to the emission reduction in gaseous precursors (Zhai et al., 2019; Yan et al., 2021).

The spatial distributions of OC and EC exhibit similarities (Fig. 10d1 and e1), consistent with the finding of a strong correlation between OC and EC in winter (Cao et al., 2007). The low temperature and weakened photochemical reactions in winter reduced secondary OC (SOC) generation, and primary OC (POC) and EC mainly originate from direct anthropogenic emissions, such as combustion (Guo, 2016). Compared to the OBS, the mass concentrations in FR (Fig. 10d2–d3) and FOR (Fig. 10e2–e3) are significantly overestimated over a wide range. Similar overestimations have also been reported in the global reanalysis datasets CAMS and MERRA-2, likely attributed to the hygroscopic growth scheme of carbonaceous aerosols in the models, poorly constrained semi-volatile species escaping from primary organic aerosols (Soni et al., 2021), and aging mechanisms in the models (Huang et al., 2013). After DA, the spatial distribution of the ANA aligns entirely with that of the OBS (Fig. 10d4 and e4), with improvements in all overestimations (Fig. 10d5 and e5), and the average biases of both OC and EC at DA sites significantly decrease to 0.14 µg m⁻³ (Fig. 11d3 and e3). The VE sites show similar results to the DA sites, with average biases of less than 2 µg m⁻³ (Fig. S5d3 and e3). Specifically, for OC (Fig. 12d1 and d2), the CORR values in FOR and ANA are 0.18–0.71 and 0.92–1.00, respectively, with RMSE values of 7.91–26.27 µg m⁻³ and 0.16–1.45 µg m⁻³, respectively. For EC (Fig. 12e1 and e2), the CORR values in FOR and ANA are 0.01–0.66 and 0.97–1.00, respectively, with RMSE values of 5.33–16.91 µg m⁻³ and 0.01–0.26 µg m⁻³, respectively. Although significant improvements are not observed in FOR at some specific DA sites, the RMSE values at all DA sites decrease by 10 %–50 % (Fig. 12d3 and e3). The CORR values of OC and EC in ANA increase by more than 30 %, with most DA sites exceeding 200 %, and the RMSE values decrease by more than 90 % (Fig. 12d4 and e4). At VE sites (Fig. S6), significant improvements in the CORR are not observed, but the RMSE values in the FOR and ANA decrease, which indicates that DA has limited benefits for whole areas but can effectively reduce biases of entire regions.

3.3 Comparison to NAQPMS-PDAF v1.0 and global reanalysis dataset

To comprehensively evaluate the competitiveness and superiority of NAQPMS-PDAF v2.0 in generating the reanalysis datasets of the PM_2.5 chemical compositions, we assimilated the mass concentrations of the five PM_2.5 chemical components from all sites (sum of DA sites and VE sites) in February 2022 to generate a reanalysis dataset. We compared our reanalysis dataset with the global reanalysis (RA) datasets (CAMSRA and MERRA-2) and NAQPMS-PDAF v1.0 output. Figure 13 illustrates the spatial distribution of the monthly average mass concentrations for the five chemical components. Compared to the OBS (Fig. 13a1 and c1), CAMSRA underestimates the NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ concentrations and fails to capture the high-value concentration in the northern part of Henan Province (Fig. 13a2 and c2). In contrast, MERRA-2 overestimates the concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC (Fig. 13b2, d2, and e2), particularly SO${}_{\mathrm{4}}^{\mathrm{2}-}$, exhibiting a large region with inaccurately high concentrations. Moreover, CAMSRA (approximately 80×80 km²) and MERRA-2 (55×70 km²) have significantly lower spatial resolutions compared to NAQPMS-PDAF v2.0 (5×5 km²). Therefore, NAQPMS-PDAF v2.0 provides a more detailed description of the pollution characteristics of chemical components in northern China and the surrounding areas compared to RA.

Figure 13 Spatial distribution of the monthly averaged concentration of five PM_2.5 chemical components for observations (OBS; a1–e1), global reanalysis (RA) data (a2–e2), NAQPMS-PDAF v1.0 analysis data (a3–e3), and NAQPMS-PDAF v2.0 analysis data (a4–e4).

Although NAQPMS-PDAF v1.0 demonstrates a superior spatial representation of the five chemical components when compared to RA, it fails to capture the high-value concentrations of NH${}_{\mathrm{4}}^{+}$ in the northwest of Henan Province and correct the high-value concentrations of NH${}_{\mathrm{4}}^{+}$ in the central and western areas of Hebei Province (Fig. 13a3). Moreover, the scattered high-value concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ in the North China Plain do not align with the spatial characteristics of the OBS (Fig. 13b3). Notably, NAQPMS-PDAF v1.0 exhibits poor performance in interpreting OC and EC, with significant overestimations in a wide range (Fig. 13d3 and e3), which indicates that NAQPMS-PDAF v1.0 is weaker than NAQPMS-PDAF v2.0 in terms of DA performance on chemical components, primarily due to insufficient propagation of observations. In NAQPMS-PDAF v2.0, the LKNETF algorithm with an adaptive forgetting factor is more suitable for the nonlinear and non-Gaussian situations compared to EnKFs in NAQPMS-PDAF v1.0, and the ensemble perturbation with a non-Gaussian distribution can better represent the reasonable error distribution of model states.

Table 3 presents a quantitative comparison of three reanalysis datasets. Compared to the CORR of NAQPMS-PDAF v2.0 (0.86–0.99), the CORR of RA for the five chemical components is significantly lower (0.42–0.55). Moreover, NAQPMS-PDAF v1.0 exhibits significantly poorer consistency in SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC, with CORR values ranging from 0.35 to 0.57. NAQPMS-PDAF v2.0 has lower overall RMSE values (0.14–3.18 µg m⁻³) compared to RA and NAQPMS-PDAF v1.0, with RMSE values ranging from 4.51 to 12.27 µg m⁻³ and 2.46 to 15.50 µg m⁻³, respectively. The characteristics of R² are similar to those of CORR and RMSE. For NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$, NAQPMS-PDAF v2.0 (0.85 and 0.93) and v1.0 (0.80 and 0.96) are much higher than RA (0.09 and 0.13). Notably, for SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC, NAQPMS-PDAF v2.0 (0.74–0.98) is significantly higher than v1.0 (−0.16–0.25) and RA (−0.15–0.25). Overall, NAQPMS-PDAF v2.0 more accurately and consistently interprets the five chemical components, particularly for NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC. The reasons are summarized as follows. (1) The DA frequency of CAMSRA is 12 h, which is lower than the hourly DA frequency in NAQPMS-PDAF v2.0. (2) CAMSRA only assimilates satellite retrievals (Inness et al., 2019), and MERRA-2 only assimilates aerosol optical depth (AOD) from both ground-based and space-based remote sensing platforms (Randles et al., 2017). The aerosol optical information analysis increment cannot be allocated to each chemical component accurately and reasonably due to the lack of a deterministic relationship between aerosol optical information and PM_2.5 chemical components. (3) NAQPMS-PDAF v1.0 has evident DA shortcomings for chemical components due to the limited DA algorithm under the assumption of a linear model or system, inappropriate ensemble perturbation under the assumption of Gaussian distribution, and inadequate observational modules. (4) The state variable structure in NAQPMS-PDAF v1.0 cannot effectively mitigate the impact of spurious correlations between chemical component variables, even when using analytical localization.

Table 3 Statistical indicators (CORR, RMSE, and R²) of five PM_2.5 chemical components for global reanalysis (RA) data, NAQPMS-PDAF v1.0 analysis data, and NAQPMS-PDAF v2.0 analysis data.

3.4 The uncertainty in NAQPMS-PDAF v2.0

In ensemble DA, the ensemble members represent possible values of the model states, and the ensemble sampling can determine the uncertainties in the model states. Therefore, the ensemble generation directly affects the propagation of observations and subsequently impacts the final DA performance. Previous studies have generated ensemble members based on the uncertainties in emission species and the Gaussian distribution assumption to satisfy the requirements of EnKF algorithms (Kong et al., 2021; Wang et al., 2022). However, the true error probability distribution of emission species is not an ideal Gaussian distribution, and the assumption will introduce errors. In this study, we coupled the hybrid nonlinear DA algorithm (LKNETF) with NAQPMS to handle the nonlinear and non-Gaussian situations, which combines the stability of LETKF with the nonlinearity of LNETF. Therefore, we evaluate the performance of ensemble members with different uncertainties and error probability distributions in NAQPMS-PDAF v2.0 through two groups of sensitivity experiments.

The first group of experiments (T1–T5) involves controlling the SO₂ uncertainty as a fixed value of 200 % and transforming the distribution of the perturbation coefficient matrix. The second group of experiments (M1–M5) focuses on assessing the influence of SO₂ uncertainty on NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ DA based on a fixed non-Gaussian distribution (m3=1, m4=6). Figure 14 shows the statistical indicators of the five chemical components under different error probability distributions, including a Gaussian distribution (T1) and four non-Gaussian distributions (T2–T5). The mean CRPS and RMSE in T2 and T4 are lower than those in T1, T3, and T5, and the CORR values in T2 and T4 are higher than those in T1, T3, and T5, indicating that the DA performance of the non-Gaussian distribution assumption is superior to that of the Gaussian distribution assumption. Moreover, positively skewed non-Gaussian distribution performs better than negatively skewed distribution. Except for SO${}_{\mathrm{4}}^{\mathrm{2}-}$, the performance in T2 outweighs that in T4 for other chemical components, implying that higher kurtosis harms the chemical component DA.

Figure 14 Statistical indicators (mean CRPS a, RMSE b, and CORR c) of five PM_2.5 chemical components for five perturbation experiments based on distribution.

SO₂ is a crucial precursor of NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$, and perturbing SO₂ affects the forecast and simulation of NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$. Table 4 presents statistical indicators of NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ analysis fields based on ensemble perturbations with different SO₂ uncertainties (12 %–300 %). Increasing the uncertainty in SO₂ from 12 % to 200 % leads to a decrease in the mean CRPS in the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ analysis field from 2.67 to 1.40, an increase in the CORR from 0.51 to 0.74, and a reduction in the RMSE from 4.10 to 2.37 µg m⁻³. Similarly, the mean CRPS in the NH${}_{\mathrm{4}}^{+}$ analysis field decreases from 0.98 to 0.77, the CORR increases from 0.88 to 0.91, and the RMSE decreases from 1.55 to 1.33 µg m⁻³. This indicates that increasing the uncertainty in SO₂ improves the DA performance on NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ because the higher SO₂ uncertainty makes SO₂ sufficiently perturbed, and the estimated error probability distribution is closer to the real distribution, resulting in a sufficient spread of observations. However, when the uncertainty in SO₂ reaches 300 %, the statistical indicators do not significantly improve and even worsen because excessively high SO₂ uncertainty causes the estimated error probability distribution to deviate from the true distribution. Thus, selecting appropriate uncertainties for emission species is crucial in aerosol chemical component DA.

To summarize, the non-Gaussian distribution assumption outperforms the Gaussian distribution assumption in NAQPMS-PDAF v2.0. Positive skewness performs better than negative skewness, and excessively high kurtosis should be avoided. Additionally, appropriately increasing the uncertainty in SO₂ enhances the DA performance on NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$. Future studies should conduct more sensitivity experiments on emission species perturbation to determine suitable schemes for different aerosol chemical components.

Table 4 Statistical indicators (mean CRPS, CORR, and RMSE; µg m⁻³) of five PM_2.5 chemical components for five perturbation experiments based on SO₂ emission uncertainty.

4 Conclusions

In this paper, we online coupled NAQPMS with PDAF-OMI to develop a novel hybrid nonlinear DA system (NAQPMS-PDAF v2.0) with level-2 parallelization based on a hybrid localized Kalman–nonlinear ensemble transform filter (LKNETF) for the first time. Compared to NAQPMS-PDAF v1.0, NAQPMS-PDAF v2.0 with OMI can be applied with multiple component types and nonlinear/non-Gaussian situations in chemical analysis to effectively interpret five PM_2.5 chemical components (NH${}_{\mathrm{4}}^{+}$, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, OC, and EC), which has not been achieved in previous studies. The background error covariance was calculated by ensemble perturbation based on adaptive uncertainties and a non-Gaussian distribution assumption of emission species. The DA experiments were conducted based on 33 observational sites in northern China and the surrounding areas.

NAQPMS-PDAF v2.0 with LKNETF can maintain high accuracy and reliability in ensemble DA with an ensemble size of 10, which is smaller than the traditional minimum of 20 ensemble members, as observed in prior ensemble assimilation studies. The free-running (FR) fields without DA have a poor consistency with the observations, with the CORR values ranging from 0.32–0.56 and the R² values being less than 0.3, showing that SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC are significantly overestimated, while NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ are underestimated. A significant improvement was observed in the analysis (ANA) fields at the DA sites. The CORR values are not less than 0.86; the RMSE and MAE values do not exceed 3.23 and 1.49 µg m⁻³, respectively; and R² is not less than 0.74. Specifically, the CORR values for NO${}_{\mathrm{3}}^{-}$, OC, and EC are not less than 0.96, and R² is not less than 0.93. The error distributions of the five chemical components concentrate to 0, with the mean bias ranging from 0 ± 0.08 to 1.02 ± 3.07 µg m⁻³. These improvements are also found in the ANA fields at VE sites, indicating an excellent DA performance of NAQPMS-PDAF v2.0.

The ability of NAQPMS-PDAF v2.0 to interpret the spatiotemporal characteristics of the five chemical components was examined. For temporal variations, compared to the FR and forecast (FOR) fields, the ANA closely aligned with the observations (OBS) and accurately captured the peak concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ in specific periods (such as 25 February), indicating good consistency and accurate characterization. Specifically, the CORR of the ANA at the six representative sites increased by 13.64 %–89.58 % and 17.19 %–75.00 %, while the RMSE decreased by 56.03 %–83.13 % and 40.74 %–72.20 %. For spatial distributions, after DA, both NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ with positive analysis increments exhibited significant improvements in CORR and RMSE, as most DA sites showed improvements of over 150 % in CORR and over 50 % in RMSE. SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC with negative analysis increments were also improved. For OC and EC in particular, the improvements in CORR and RMSE at most DA sites were over 200 % and over 90 %, respectively. The improvements at VE sites were also identified. Consequently, DA successfully aligned the spatiotemporal characteristics of the ANA with OBS and significantly reduced the biases of five chemical components.

Compared to the global reanalysis datasets (CORR: 0.42–0.55, RMSE: 4.51–12.27 µg m⁻³) and NAQPMS-PDAF v1.0 (CORR: 0.35–0.98, RMSE: 2.46–15.50 µg m⁻³), NAQPMS-PDAF v2.0 (CORR: 0.86–0.99, RMSE: 0.14–3.18 µg m⁻³) has significant superiority in generating the reanalysis datasets of the PM_2.5 chemical compositions with high spatiotemporal resolution. Moreover, NAQPMS-PDAF v1.0 cannot capture the high-value concentrations and exhibits poor performance when interpreting SO${}_{\mathrm{4}}^{\mathrm{2}-}$, OC, and EC with CORR values ranging from 0.35 to 0.57. In contrast, NAQPMS-PDAF v2.0 interprets the five chemical components more accurately and consistently.

Finally, the uncertainties in NAQPMS-PDAF v2.0 are examined by identifying the influence of ensemble generation on ensemble DA performance. The non-Gaussian distribution assumption outperforms the Gaussian distribution assumption in NAQPMS-PDAF v2.0. Positive skewness performs better than negative skewness, and excessively high kurtosis should be avoided. Additionally, appropriately increasing the uncertainty in SO₂ enhances the DA performance on NH${}_{\mathrm{4}}^{+}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$. Future studies should conduct more sensitivity experiments on emission species perturbation to determine the schemes that are suitable for different aerosol chemical components.

The novel hybrid nonlinear DA system (NAQPMS-PDAF v2.0) can be effectively applied in the interpretation of chemical components and outperforms the reanalysis datasets in generating the five PM_2.5 chemical components with high accuracy and high consistency, thus providing a sufficient channel to investigate spatiotemporal characteristics, identify regional transport, and prevent and control aerosol composition pollution. In future work, we plan to research the vertical DA of chemical components, introduce more vertical information from more observational platforms, and verify the simultaneous DA performance of surface and vertical mass concentrations.