Operational chemical weather forecasting with the ECCC online Regional Air Quality Deterministic Prediction System version 023 (RAQDPS023) – Part 1: system description

The online version of the Regional Air Quality Deterministic Prediction System (RAQDPS) is a chemical weather forecast system that has been employed operationally by Environment and Climate Change Canada (ECCC) since 2009. It is run twice per day to produce 72 h forecasts of hourly 10 km abundance fields of three key predictands, NO₂, O₃, and PM_2.5 total mass, as well as other gas-phase chemical species, PM_2.5 chemical components, and dry and wet deposition for Canada, the contiguous U.S., and northern Mexico. The forecasts of NO₂, O₃, and PM_2.5 are needed to calculate the Air Quality Health Index (AQHI), which is used to communicate current and forecasted pollutant levels to the Canadian public. Version 023 of the RAQDPS (RAQDPS023) went into service at ECCC in December 2021 and was replaced by the RAQDPS025 in June 2024. This paper provides the first full description of any version of the online RAQDPS. After giving a brief history of the ECCC operational air quality forecasting program, we provide a comprehensive description of the RAQDPS023 forecast system as well as shorter descriptions of several upstream and downstream forecast and analysis systems. The latter include two upstream operational meteorological forecast systems that were based on version 5.1.0 of the ECCC Global Environmental Multiscale (GEM) numerical weather prediction model, one which used a global configuration, the Global Deterministic Prediction System (GDPS 8.0.0), and the other which used a regional configuration, the Regional Deterministic Prediction System (RDPS 8.0.0). An emissions processing system, an Updateable Model Output Statistics-based system for bias-corrected station-specific pollutant concentration forecasts (UMOS-AQ), and a regional objective analysis system for surface pollutant concentration fields, the Regional Deterministic Air Quality Analysis system (RDAQA 2.0.0), are also described.

The RAQDPS023 itself consisted of version 3.1.0.0 of the GEM-Modelling Air quality and CHemistry (GEM-MACH) chemistry module, which was embedded with one-way coupling within GEM 5.1.0, its meteorological host model. The meteorological configuration of the RAQDPS023 closely followed that of the RDPS 8.0.0. Details covered in this paper include a summary of the dynamical representations and physical parameterizations used in the three GEM-based forecast systems, which are highly harmonized, the chemical species and parameterizations used in the MACH chemistry module, including gas-phase, aqueous-phase, and inorganic heterogeneous schemes and associated numerical solvers, system inputs, including both anthropogenic and natural emissions of chemical species, system outputs, and run configuration, strategies, and timings. One simplification in addition to the use of the condensed ADOM-2 gas-phase chemistry scheme that was made to reduce RAQDPS023 execution time for operational deployment was to represent the particulate matter (PM) size distribution with only two aerosol particle size bins, one corresponding to particle diameters in the 0–2.5 µm range (“fine particles” or PM_2.5) and the other to the 2.5–10 µm range (“coarse fraction” or PM_cf). A second simplification was to represent the chemical composition of PM_2.5 with only nine chemical components, and a third simplification was to use a longer time step (900 s) for the time integration of atmospheric chemistry than the time step used for time integration of atmospheric dynamics and physics (300 s). Even so, activating the MACH module increased RAQDPS023 run time by a factor of 4.4 on average compared to meteorology only, partly due to the cost of the integration of chemistry but partly to the increased cost of the GEM dynamical core due to the advection with imposed shape preservation and mass conservation of 57 additional chemical tracers. The role of the RAQDPS-FW023, a second chemical weather forecast system that was identical to the operational RAQDPS023 (or RAQDPS-OP023) except for the addition of near-real-time biomass burning emissions, is also described. Biomass burning emissions for Canada and the U.S. estimated from satellite measurements were first calculated by the Canadian Forest Fire Emissions Prediction System (CFFEPS) version 4.1 before each RAQDPS-FW023 run was launched. Outputs from the two RAQDPS versions were then used to produce forecasts of wildfire smoke transport and diffusion. The paper closes by summarizing the key upgrades made to the RAQDPS025, the current version of the ECCC operational chemical weather forecast system, and then describing some possible future improvements and updates. A companion paper by Moran et al. (2026) presents the results of a comprehensive, five-year performance evaluation of prospective and retrospective annual air quality simulations made with the RAQDPS023.

1 Introduction

Forecasting tomorrow's weather means being able to predict the short-term future dynamical and physical state of the atmosphere. Operational weather forecasting, also known as numerical weather prediction (NWP), began in the 1950s once electronic computers started to become available commercially (e.g., Lorenz, 2006; Bauer et al., 2015). The chemical state of the atmosphere, on the other hand, particularly the concentrations of air pollutants of concern for human health, is often referred to as air quality (AQ). Recent epidemiological studies have determined that air pollution is the fourth leading cause of early death globally (Boogaard et al., 2019; Murray et al., 2020; WHO, 2021). However, predicting tomorrow's air quality, or chemical weather as it is also known, is even more difficult than predicting tomorrow's weather. As a consequence, operational regional AQ forecasting did not become possible until about 2000 (McHenry et al., 2004; Kukkonen et al., 2012; Zhang et al., 2012a), roughly a half century after the start of operational weather forecasting.

Air quality was initially only a local concern in the vicinity of large sources of air pollutants, including major urban centres (Pasquill and Smith, 1983). In the 1970s, however, emerging recognition of human health impacts caused by regional photochemical smog and the environmental harm caused by regional acidic deposition (or “acid rain”) spurred many laboratory and field studies to better understand these larger-scale AQ issues (e.g., Hecht et al., 1974; Likens and Bormann, 1974; Eliassen and Saltbones, 1975; Barrie and Georgii, 1976; Calvert et al., 1978; Harrison et al., 1978; Schiermeir, 1978). The development of the first prognostic, three-dimensional, regional AQ models capable of simulating photochemical smog, especially ozone, and acidic deposition then followed in the 1980s, building on the steady increase in our scientific understanding of these issues (e.g., Chang et al., 1987; Seinfeld, 1988; Venkatram et al., 1988; Russell and Dennis, 2000). At first these regional AQ models, which are also known as chemical transport models (CTMs) and atmospheric chemistry models (Kukkonen et al., 2012), were only applied to retrospective case studies for research and policy applications (e.g., Clark et al., 1989; Fung et al., 1991; Roselle and Schere, 1995; Dunker et al., 1996; Middleton, 1997). Routine regional AQ forecasting was held back by two factors – the order-of-magnitude larger computational burden of AQ models relative to NWP models and the run-time constraints imposed by the need to disseminate forecasts in a timely fashion (McHenry et al., 2004) – but continuing improvements in computer capabilities have facilitated the introduction and rapid expansion of the use of operational AQ forecast systems in North America, Europe, and Asia over the past two decades (e.g., Grell and Baklanov, 2011; Kukkonen et al., 2012; Zhang et al., 2012a, b; Brasseur et al., 2019; WMO, 2020; Brasseur and Kumar, 2021).

The majority of AQ forecast systems employ a so-called “offline” model architecture, where an NWP model is run first and then a CTM is run using forecast meteorological fields from the NWP model as inputs. More recently, however, some AQ forecast systems have adopted an “online” architecture where meteorological and chemical variables are predicted simultaneously, in effect a merged NWP-CTM system (e.g., Grell et al., 2005; Zhang, 2008; Grell and Baklanov, 2011; Baklanov and Zhang, 2020).

Canada was one of the first countries to implement an operational regional AQ forecast system (e.g., Zhang et al., 2012a). Environment and Climate Change Canada (ECCC), Canada's federal environment ministry, has run an operational, continental-scale Regional Air Quality Deterministic Prediction System (RAQDPS) for North America since 2001. The earliest version of the RAQDPS employed an offline CTM called CHRONOS (Canadian Hemispheric Regional Ozone and NO_x System). This system was run once per day on a 21 km horizontal grid over North America with 25 vertical levels from the surface to 6 km to make 48 h forecasts of surface ozone (O₃) volume mixing ratio (Pudykiewicz et al., 1997; Pudykiewicz and Koziol, 2001; McHenry et al., 2004; McKeen et al., 2005; Ménard and Robichaud, 2005; Moran et al., 2013; Robichaud and Ménard, 2014). The meteorological “driver” (i.e., NWP model) used with CHRONOS was the limited-area ECCC Regional Deterministic Prediction System (RDPS), the operational, nested, regional configuration of the ECCC Global Environmental Multi-scale (GEM) NWP model (Côté et al., 1998a, b; Fillion et al., 2010; Caron et al., 2015). The surface concentration of particulate matter (PM) with aerodynamic diameter smaller than 2.5 µm (PM_2.5) was added as a predictand to CHRONOS in 2003 (Pudykiewicz et al., 2003; McKeen et al., 2007).

In 2009 the offline CHRONOS CTM in the RAQDPS was replaced by the online, one-way-coupled GEM-MACH (GEM–Modeling Air quality and CHemistry) chemical weather model, which included improved representations of gas-phase, aqueous-phase, and aerosol chemistry in the new MACH chemistry module that was added to the GEM model. This new generation of the RAQDPS used a new, continental-scale 15 km horizontal grid with 58 vertical levels extending from the surface to 0.1 hPa (Anselmo et al., 2010; Moran et al., 2010, 2011). One design goal for the development of this next generation was to align its meteorological component and configuration with the operational RDPS as closely as possible to ensure equivalent meteorological forecasts from these two closely related systems.

From 2009 to 2025 there have been 24 upgrades of varying magnitude made to the GEM-MACH-based version of the RAQDPS to improve forecast skill and to adjust to new versions of the operational RDPS weather forecast model and upgrades to ECCC's computer system hardware and software. For example, version 023 of the RAQDPS (RAQDPS023), the subject of this paper, employed a continental forecast grid with 10 km horizontal grid spacing and 84 hybrid vertical levels from the surface to 0.1 hPa. As a second example, a key input to any AQ forecast model is a detailed representation of current anthropogenic and natural emissions of air pollutants and their precursors: the input emissions files used by the RAQDPS have been updated eight times since 2009. Other significant changes have included numerous improvements and some bug fixes to the chemistry module code, several modifications to the vertical discretization, the introduction of a duplicate version of the RAQDPS that considers near-real-time (NRT) wildfire emissions, and an extension of the forecast period from two days to three days. Table A1 lists all of these upgrades and indicates the biggest ones, four of which involved major updates to the input anthropogenic emissions files resulting from the availability of newer national emission inventories.

The operational RAQDPS023, also referred to as the RAQDPS-OP023, went into service on 1 December 2021 (Moran et al., 2021b). It was run routinely twice daily until 10 June 2024 to produce 72 h forecasts of hourly surface abundance fields of O₃, nitrogen dioxide (NO₂), PM_2.5, and other chemical species over North America. These first three species were needed to calculate Canada's national Air Quality Health Index (AQHI), a health-based, multi-pollutant, additive, no-threshold AQ index (Sect. 5.1). The RAQDPS-FW023, a duplicate version of the RAQDPS-OP023 that also includes NRT wildfire emissions but is otherwise identical, also went into service on 1 December 2021 (Chen and Menelaou, 2021). Forecasts made by this second system allowed the location and evolution of wildfire plumes to be calculated as the difference between forecasts from the two model versions (Pavlovic et al., 2016; Munoz-Alpizar et al., 2017; Chen et al., 2019b).

The goal of this paper is to describe and document all components of the RAQDPS-OP023 and RAQDPS-FW023 forecast systems. As such it is the first publication to provide a comprehensive and detailed description of any version of the online, GEM-MACH-based RAQDPS. A companion paper (Moran et al., 2026) presents a detailed performance evaluation of five years of RAQDPS023 simulations consisting of the first year of RAQDPS-OP023 forecasts and four years of retrospective annual simulations for the 2013–2016 period that used year-specific emissions. It should be noted that the RAQDPS023 was replaced in June 2022 by the RAQDPS024, a new but equivalent version that was necessary due to the migration of all ECCC forecast systems to a new computer system. The operational GEM-MACH code and the input emissions were unchanged between these two versions. Another release, the RAQDPS025, was implemented operationally on 11 June 2024 (see Table A1), but despite a number of upgrades introduced for the RAQDPS025 that are described near the end of this paper, almost all of the details related to the RAQDPS023 that are presented in this paper also apply to the RAQDPS025.

The rest of the paper is organized as follows. Section 2 provides a brief overview of the meteorological aspects of the RAQDPS023 operational chemical weather forecast system and the two upstream operational NWP forecast systems that provide the RAQDPS023 with meteorological initial conditions and boundary conditions. These two upstream systems are global and regional configurations of the GEM NWP model. Section 3 provides a detailed description of the MACH component of the RAQDPS023, including its chemical parameterizations and input anthropogenic and modelled natural emissions. Section 4 describes computational aspects of the RAQDPS023, including time integration and process splitting, chemical initial conditions and boundary conditions, advection scheme and tracer shape preservation and mass conservation, computer hardware and code parallelization, and forecast run strategies. Section 5 describes RAQDPS023 outputs, post-processing and downstream systems, and routine operational performance evaluation. Finally, Sect. 6 describes updates made for the current RAQDPS025 operational forecast system and some possible future improvements, and Sect. 7 summarizes the paper.

2 Meteorological forecast systems

The meteorological forecast component of the RAQDPS023 is the GEM NWP model (Côté et al., 1998a, b; Girard et al., 2014; McTaggart-Cowan et al., 2019a; Ritchie et al., 2022). As the GEM model has been described in more detail in many other publications, only its main characteristics are summarized here with an emphasis on those parameterizations that are also important for chemical weather forecasting. These include the treatment of vertical fluxes from land and water surfaces, turbulent mixing in the surface layer and the planetary boundary layer (PBL), shortwave and longwave radiation, and clouds and precipitation.

The GEM model can be configured for either a global domain or a regional domain (Côté et al., 1998a, b). The ECCC operational medium-range global forecast system based on GEM is named the Global Deterministic Prediction System (GDPS) while the operational short-range regional forecast system for the North American region is named the Regional Deterministic Prediction System (RDPS) (Buehner et al., 2015; Caron et al., 2015; McTaggart-Cowan et al., 2019a; Ritchie et al., 2022). The RDPS is also referred to as a limited-area model (LAM). These two meteorological forecast systems are very similar in terms of their dynamical cores (see Sect. S1.1 of the Supplement to this paper), but additional efforts were also made in recent years to harmonize them further in terms of numerics and physical parameterizations both for consistency and to reduce development and maintenance burdens. For example, both the GDPS and RDPS now employ the same number of vertical levels, the same model lid height, and the same set of atmospheric physics parameterizations (Mailhot et al., 1998; McTaggart-Cowan et al., 2019a).

New operational versions of the GDPS and RDPS, GDPS 8.0.0 and RDPS 8.0.0, were introduced on 1 December 2021 at the same time as the RAQDPS023 (CMC-GDPS-8.0.0, 2021a, b; CMC-RDPS-8.0.0, 2021a, b; Moran et al., 2021b; CMC-RAQDPS-023, 2021). Given that the focus of this paper is a regional chemical weather forecast system, the main reason that the GDPS 8.0.0 global system was included here is that the RDPS 8.0.0 regional system was one-way coupled with GDPS 8.0.0 forecasts so that the regional NWP system was dependent on the global NWP system. In addition, the meteorological aspects of the RAQDPS023 chemical weather forecast system were identical to the RDPS 8.0.0 configuration, and the RAQDPS023 was one-way coupled to the RDPS8.0.0 for meteorology since the RAQDPS023 horizontal domain was embedded within the RDPS8.0.0 horizontal domain (Sect. 4.2). Figure 1 shows the data-flow relationships between the GDPS 8.0.0, the RDPS 8.0.0, and the RAQDPS023, namely a sequence of three linked operational forecast systems where the GDPS 8.0.0 and RDPS 8.0.0 are “upstream” dependencies of the RAQDPS023.

Figure 1 RAQDPS023 “upstream” and “downstream” dependencies on other ECCC operational systems.

All three systems were based on GEM version 5.1.0. Most of the following description of GEM 5.1.0 thus applies to all three systems, but differences are noted where they occur. The main differences between the global and the two regional GEM-based systems were related to run configurations, namely domain size, horizontal grid spacing, integration time step, initialization, and boundary conditions.

2.1 Grids, coordinate systems, and nesting

The GDPS 8.0.0 employed an overset Yin–Yang global horizontal grid system that combines two perpendicular and overlapping, rotated latitude-longitude limited-area horizontal grids referred to as the Yin grid and the Yang grid. The globe is thus divided into two regions that look somewhat reminiscent of the covering and seams of a baseball (Qaddouri and Lee, 2011). This global grid system has two major benefits: it avoids pole-related singularities and convergence issues and it is suitable for use with massively parallel computers (Kageyama and Sato, 2004). Forecasts were made independently at each time step on the two regional LAM grids and were then reconciled across the overlap region (Qaddouri and Lee, 2011). The horizontal grid spacing for both the Yin and Yang grids used by the GDPS 8.0.0 was quasi-uniform 0.135° (∼ 15 km).

The rotated latitude-longitude LAM horizontal grid of the RDPS 8.0.0 covered the North American continent and adjacent oceans but was entirely a subdomain of the GDPS's Yin grid (Fig. 2). Horizontal grid spacing was quasi-uniform 0.09° (∼ 10 km). The RAQDPS023 horizontal grid, also shown in Fig. 2, was an embedded subgrid of the RDPS 8.0.0 grid with 0.09° grid spacing and grid-point superposition (i.e., co-location). The RDPS 8.0.0 horizontal grid was 1108 × 1082 in size vs. 772 × 642 for the RAQDPS023 horizontal grid. The smaller domain of the chemical weather model was necessary to balance the additional computational burden of the MACH chemistry module (see Sect. 4.1).

Figure 2 Horizontal domains of RAQDPS023 (blue inner area) and RDPS 8.0.0 (surrounding red area plus blue inner area).

In terms of finite differencing, the horizontal discretization used for all three model grids was an Arakawa C grid, which Côté et al. (1998a) argued is more suitable for massively parallel computer architectures and for mesoscale applications. The vertical discretization for all three model grids was a staggered Charney–Phillips grid that was used with a log-hydrostatic-pressure-type ζ (or hybrid) vertical coordinate (Holdaway et al., 2013a, b; Girard et al., 2014). The GDPS 8.0.0, RDPS 8.0.0, and RAQDPS023 also all used the same 84 hybrid vertical levels stretching from the Earth's surface to 0.1 hPa. The lowest momentum “full” hybrid levels were located at approximately 20, 61, 114, 181, and 263 m a.g.l., and the lowest thermodynamic “half” hybrid levels were located at approximately 10, 40, 87, 147, and 222 m a.g.l. Vertical layer thickness was monotonic increasing with distance from the Earth's surface to the upper troposphere (where it then decreased to better resolve the tropopause region), and there were 11 thermodynamic hybrid levels located within the first kilometer above the surface.

GDPS 8.0.0 forecasts started before RDPS 8.0.0 forecasts so that one-way nesting could be used to supply meteorological lateral boundary conditions (LBCs) from the global weather forecasts to the LAM regional weather forecasts (Benoit et al., 1997; Ritchie et al., 2022; see also Sect. 4.2). Similarly, RDPS 8.0.0 forecasts started before RAQDPS023 forecasts so that one-way nesting could again be used to supply meteorological LBCs from the RDPS 8.0.0 LAM forecasts to the RAQDPS023 LAM forecasts (see Figs. 1 and 2 and Sect. 4.2). This one-way nesting technique is sometimes referred to as “piloting”. In fact, successive versions of the RAQDPS have been “piloted” or “driven” by successive versions of the RDPS since 2009 (Table A1).

2.2 RAQDPS023 meteorological configuration

Table 1 provides a concise summary of the meteorological components and attributes of the RAQDPS023, including the dynamical core and various physical parameterizations. References are also provided for each component. These components and attributes were identical to those used by the RDPS 8.0.0 with two exceptions: horizontal domain size and LBC source (cf. CMC-RDPS-8.0.0, 2021b). For further information about Table 1, McTaggart-Cowan et al. (2019a) provide an excellent recent summary of ECCC implementations of the physical parameterizations listed in Table 1, Sect. S1 also provides some information about many of these components and attributes, and Sect. 4 covers some computational aspects.

Table 1 RAQDPS023 meteorological configuration and references for GEM host model (see also CMC-RAQDPS-023, 2021).

3 Chemical forecast module

The MACH (Modelling Air quality and CHemistry) module is an online chemistry module embedded in the GEM NWP model code to produce the GEM-MACH chemical weather model. It predicts the evolution in time of concentrations and removal rates of a number of gas-phase chemical species, including O₃ and NO₂, as well as size- and composition-resolved PM. Many of the chemistry parameterizations used by the MACH module were adapted from two earlier ECCC atmospheric chemistry codes, the Canadian Aerosol Module (CAM), which was designed to represent aerosol processes in climate models (e.g., Gong et al., 1997a, b, 2002, 2003), and A Unified Regional Air Quality Modelling System (AURAMS) CTM, which was an offline AQ model with a size- and composition-resolved representation of PM developed by ECCC for research and policy applications (e.g., Moran et al., 1998; Zhang et al., 2002a; Gong et al., 2006; Stroud et al., 2008; Makar et al., 2009; Levy et al., 2010). Furthermore, as noted below several chemistry parameterizations used by AURAMS were adapted from an earlier regional CTM used by ECCC, the Acid Deposition and Oxidant Model (ADOM: Venkatram et al., 1988, 1992; Misra et al., 1989; Fung et al., 1991, Karamchandani and Venkatram, 1992; Macdonald et al., 1993; Li et al., 1994).

3.1 Alignment with weather forecast models

For the RAQDPS023 the GEM-MACH code version was updated from version 3.0.2 to version 3.1.0.0 (see Table A1 and Moran et al., 2021b). This updated version incorporates the GEM v5.1.0 code that was also used by the GDPS 8.0.0 and RDPS 8.0.0 (CMC-GDPS-8.0.0, 2021a; CMC-RDPS-8.0.0, 2021a).

As noted in Sect. 2.1 the RAQDPS023 meteorological configuration was also very closely aligned with that of the RDPS 8.0.0 regional weather forecast model. First, the RAQDPS023 employed (i) the same set of dynamics options and physics parameterizations as the RDPS 8.0.0, including those for advection, vertical turbulent diffusion, radiation, grid-scale and subgrid-scale (SGS) clouds and precipitation, and gravity wave drag and orographic blocking (Table 1), (ii) the same vertical grid and a closely related horizontal grid (the RDPS horizontal grid domain is larger; see Sect. 2.1 and Fig. 2), and (iii) the same dynamics/physics integration time step (Table 1). As a consequence, weather forecasts made by the two systems are nearly identical over the RAQDPS horizontal domain.

Second, the chemistry parameterizations used by the RAQDPS023 made use of many meteorological fields predicted by GEM. For example, the treatment of non-chemical cloud processes resides in the GEM physics module, which provides a number of meteorological fields to the MACH chemistry module, including cloud water mixing ratio (liquid and solid), precipitation production rate (auto-conversion and coalescence), precipitation evaporation rate, and precipitation vertical fluxes (both liquid and solid) (Mailhot et al., 1998; Gong et al., 2006, 2015; McTaggart-Cowan et al., 2019a). Table 2 summarizes the meteorological fields predicted by GEM that are needed by various MACH parameterizations in the RAQDPS023.

Table 2 List of GEM-predicted meteorological fields needed by the MACH module. Two-dimensional fields are associated with the Earth's surface while three-dimensional fields are associated with the atmospheric column.

Third, as noted in Table 1 the GEM dynamical core used by the RAQDPS023 included two meteorological tracers: water vapour and cloud condensate. MACH considers a much larger number of chemical tracers, but these additional tracers are implemented in the GEM-MACH code in the same way as the meteorological tracers. In fact, MACH chemistry is essentially a run-time option of GEM that can either be activated to run GEM-MACH for chemical weather forecasts or not activated so that just GEM is run to make meteorological forecasts. Moreover, because GEM-MACH is an online system, it makes use of GEM schemes and Fortran2003 computer code for those dynamical and physical processes that are common to both meteorological and chemical tracers, including advection and vertical diffusion. One exception, however, is SGS convection, where SGS vertical transport of chemical tracers mirroring that of meteorological tracers has not yet been implemented.

Fourth, the same a posteriori treatments for shape preservation and mass conservation applied in the RDPS 8.0.0 for the advection of meteorological tracers (Table 1) are also applied in the RAQDPS023 for the advection of both meteorological and chemical tracers. Shape preservation and mass conservation are very relevant for a regional AQ model since there can be many sharp chemical concentration gradients in the vicinity of strong local emissions sources, and the resulting loss of shape preservation after advection can introduce negative concentrations (e.g., Flemming et al., 2015).

And fifth, in order to use the treatment of vertical turbulent diffusion from GEM for chemical tracers, area emissions and dry deposition of gases and particles are treated as flux boundary conditions at the bottom boundary. Note that two small differences in the parameterization of vertical diffusion for meteorological tracers vs. chemical tracers used by the RAQDPS023 were introduced in 2019 when the RAQDPS021 became operational (Moran and Ménard, 2019). In the RAQDPS021 the height of the lowest model half-level (i.e., first thermodynamic level) was reduced from 20 to 10 m and the scheme used to diagnose PBL height was also changed to a Ri_b-based scheme (Sect. S1.4), which resulted in PBL heights as small as 2 m being predicted under some stable conditions (Moran and Ménard, 2019). To avoid the occurrence of unrealistically high surface concentrations, two ad hoc modifications were introduced. One modification was to impose a minimum PBL height value of 100 m (other groups have imposed similar minimum values: e.g., Li and Rappenglueck, 2018). The second modification was to inject surface emissions equally into the two lowest model layers rather than just in the lowest layer as before 2019, which was equivalent to maintaining the 40 m lowest full-layer thickness considered before the RAQDPS021. Both modifications were carried over to the RAQDPS023.

3.2 Representation of particulate matter

The term “particulate matter” (or aerosol particles) refers to the mixture of solid particles and liquid droplets found suspended in air. The size of these particles can range from the “ultrafine” (aerodynamic diameters less than 0.1 µm) to the “giant” (aerodynamic diameters greater than 10 µm) (e.g., Wang et al., 2014a). PM is of interest for many reasons: for example, aerosol particles can interact with radiation via scattering and absorption (aerosol direct effect); they can act as condensation and ice nuclei in clouds and interact with cloud microphysics (aerosol indirect effect); and they can enter the human body via respiration, where they may impact human health. In fact, so-called “fine” particles (PM_2.5) are a major cause of human mortality and morbidity (e.g., Lelieveld and Pöschl, 2017; Murray et al., 2020). As a consequence, many countries have legislated ambient air quality standards for PM_2.5, including Canada (https://ccme.ca/en/air-quality-report, last access: 18 March 2026) and the U.S. (https://www.epa.gov/criteria-air-pollutants/naaqs-table, last access: 18 March 2026), and the World Health Organization has issued health-based guidelines for PM_2.5 and PM₁₀ ambient concentrations (WHO, 2021).

PM is a complex pollutant: (i) it is typically composed of particles with a wide size range spanning orders of magnitude; (ii) it can contain many different elements and compounds, both inorganic and organic; (iii) it is subject to many different physical and chemical processes (see below); and (iv) it has many different types of sources, including both particle-phase primary sources and gas-phase secondary sources (e.g., Fuzzi et al., 2015). Representing size- and composition-resolved PM in an AQ model with an affordable number of tracer species is thus a significant challenge both scientifically and computationally.

Two conceptual approaches have been used in AQ models to represent the PM size distribution: a sectional representation, in which the PM size distribution is divided in n contiguous, nonoverlapping, discrete size sections or size bins (e.g., Gelbard and Seinfeld, 1980; Seigneur et al., 1986; Jacobson, 1997; Zhang et al., 1999); and a modal representation, in which the PM size distribution is assumed to be composed of n distinct populations of particles, each of whose size distribution can be described by an analytical modal distribution function (e.g., Binkowski and Shankar, 1995; Whitby and McMurry, 1997). Many AQ models with a size-resolved treatment of PM have employed a sectional representation, (e.g., Jacobson, 1997; Kukkonen et al., 2012; Zhang et al., 2012a, b; WMO, 2020). For each model, though, the number of sections (n) that are considered must be chosen.

For GEM-MACH two different values of n have been used, depending upon the application. For operational AQ forecasting with the RAQDPS, only two size bins have been used in order to reduce model execution times: PM_2.5 or “fine PM”, which corresponds to particle diameters in the 0–2.5 µm range; and PM_cf or “coarse fraction”, which corresponds to particle diameters in the 2.5–10 µm range. Taken together these two size bins constitute PM₁₀ (“coarse PM”). These two size bins were chosen to enable prediction of PM_2.5 and PM₁₀, the two PM size ranges of concern for AQ standards in North America, following the approach of an earlier size-resolved PM model (Middleton, 1997). These two size bins are also consistent with measurement studies, which have found that atmospheric aerosol volume and mass size distributions generally have only two distinct modes, an accumulation mode (0.1 to 2 µm diameter range) and a coarse mode (diameter > 2 µm), except for very clean conditions or close to sources of hot gases where a distinct (and smaller) nucleation or Aitken mode is evident (e.g., Whitby, 1978; Morawska et al., 1999). Whitby (1978) also argued that there is little interaction between particles in the accumulation and coarse modes. Note too that PM_2.5 and PM₁₀ are the two PM size ranges that are considered by the Canadian, U.S., and Mexican national emissions inventories (see Sect. 3.11.1) and by a number of North American surface measurement networks (e.g., Brook et al., 1999b; NSTC, 2013; Moran et al., 2026). However, as discussed below in Sect. 3.3.4 and 3.8 it was necessary to introduce an additional subdivision of the PM_2.5 and PM_cf size bins in order to reduce errors in the numerical solution of several aerosol process parameterizations. It is should be noted that some model applications do need the PM size distribution to be resolved in more detail, such as calculations of atmospheric visibility and aerosol optical depth or studies of meteorology-aerosol interactions through aerosol direct and indirect effects. For such applications GEM-MACH is typically run with 12 size bins (e.g., Gong et al., 2015; Makar et al., 2015a, b; Ghahreman et al., 2021; Majdzadeh et al., 2022).

Another complication related to the PM size distribution is that there are multiple definitions of aerosol particle diameter. The most common definition is aerodynamic diameter, which is defined as the diameter of a spherical particle with a density of 1 g cm⁻³ that behaves similarly to the particle of interest. This is the definition used for PM_2.5 and PM₁₀ for AQ standards, for reporting most ambient measurements, and for preparing emissions inventories. However, CTMs like GEM-MACH consider Stokes diameter, which is the diameter of a spherical particle that behaves similarly to the particle of interest and has the same density. Stokes diameter is equal to the aerodynamic diameter divided by the square root of the particle density (e.g., Seigneur and Moran, 2004), and this difference should be considered when comparing model-predicted PM values with observed PM values (see companion paper by Moran et al., 2026).

While it is important for the RAQDPS to predict PM_2.5 “bulk” mass or total mass since PM_2.5 total mass is one of the three species considered by the AQHI, it is also important for the RAQDPS to resolve the chemical composition of PM_2.5 and PM₁₀ since many PM properties and processes depend on a particle's chemical composition. In order to achieve this representation with only a modest number of chemical components so as to reduce computational costs, just nine components were considered by the RAQDPS023: sulfate (SU); nitrate (NI); ammonium (AM); elemental carbon (EC); primary organic matter (POM); secondary organic matter (SOM); crustal material (CM); sea salt (SS); and aerosol water (WA). This set of chemical components is consistent with the approach of three North American PM_2.5 speciation measurement networks (CSN, IMPROVE, NAPS) to report PM_2.5 chemical composition (see Dabek-Zlotorzynska et al., 2011; Chow et al., 2015) and with other AQ models that resolve PM composition (e.g., Table 6 of Kukkonen et al., 2012). Note, however, that all of these chemical components must be considered to be “lumped” species. In the case of particle SU, NI, and AM, they do not account for the molecular or stoichiometric form of these three inorganic constituents, which controls particle acidity. For example, particle SU may be present as sulfuric acid (H₂SO₄(aq)), ammonium sulfate ((NH₄)₂SO₄(s)), ammonium bisulfate ((NH₄HSO₄(s)), or letovicite ((NH₄)₂H(SO₄)₂(s)) (e.g., Saxena et al., 1983; Makar et al., 2003b). Particle EC, which is often referred to as black carbon (BC), is implicitly defined by the source-testing and ambient measurement methods used to measure it (e.g., Chow, 1995; Chow et al., 2015, 2018). Particle organic matter (POM and SOM) is composed of other elements such as hydrogen, oxygen, sulfur, and nitrogen, in addition to carbon (e.g., Turpin and Lim, 2001; Pang et al., 2006; Malm and Hand, 2007; Philip et al., 2014). Particle CM (or soil dust) is composed of oxides of silicon, various metals (e.g., iron, aluminum, titanium), base cations (sodium, calcium, magnesium, and potassium), and carbonate (e.g., Malm et al., 2004; Dabek-Zlotorzynska et al., 2011; Hand et al., 2017). And particle SS is a complex mixture of multiple salts and organics in addition to NaCl (e.g., White, 2008; Ault et al., 2013).

One additional aspect of PM composition is its mixing state, which refers to the variation of PM chemical composition with particle size (e.g., Heintzenberg, 1989; Jacobson, 2001; Zhu et al., 2015). At one extreme, referred to as an internal mixture, all aerosol particles of a given size at a given point in space and time are assumed to have the same chemical composition. At the other extreme, referred to as an external mixture, each aerosol particle of the same size may be composed of any one of many pure chemical species. And between these two extremes lies a very wide range of transitional mixing states. Mixing state affects a particle's hygroscopicity, organic absorptivity, and radiative properties (Stevens et al., 2022). MACH assumes that each PM size bin is internally mixed (Gong et al., 2003; Park et al., 2011; Stevens et al., 2022). This assumption has the advantage of minimizing the number of PM chemical components that must be considered since no special considerations must be made for mixing state. It is least good close to emission sources of primary PM but becomes more realistic as aerosol particles age and undergo condensation and coagulation processes (e.g., Winkler, 1973).

Size- and composition-resolved atmospheric PM was thus represented in the RAQDPS023 with two internally-mixed size bins and nine chemical components for a total of 18 size bin-chemical component tracers (e.g., PM_2.5-EC, PM_cf-NI). Sixteen of these PM tracers were prognostic fields and were advected while the two size bin-aerosol water tracers were diagnostic fields that were estimated by the Hänel (1976) scheme or the HETV code (Sect. 3.3.5 and 3.6). Although each size bin was assumed to be internally mixed, the chemical composition of the two size bins could be different. Mass mixing ratio (MMR) with units of µg size bin-component per kg dry air was used in the PM conservation-of-mass equations in MACH to describe aerosol particle abundances. Mixing ratios have the advantage of being independent of pressure and temperature. While molar (or volume) mixing ratio is an official SI unit for species abundance in air (Schwartz and Warneck, 1995), MMR is a more appropriate choice when the chemical composition of a constituent is not known or is not well-defined as is the case for “lumped” chemical species such as the above nine PM chemical components. MMR units were also used for gaseous species in the RAQDPS023 for consistency, although unit conversions may be performed before and after some process operators (e.g., gas-phase chemistry).

3.3 Aerosol microphysical schemes

A number of physico-chemical processes affect the PM population and its size distribution and chemical composition: (i) emission of primary particles of different sizes and composition; (ii) nucleation of new ultrafine particles from gas-phase precursors; (iii) growth of existing particles via condensation of non-water vapours and the shrinkage of particles via volatilization of particulate species to the gas phase (i.e., gas-particle partitioning); (iv) collision and coagulation of particles; (v) swelling and activation of particles due to water-vapour sorption processes; and (vi) size-dependent removal of particles via dry and wet deposition processes. CTMs that resolve the PM size distribution should in principle simulate all of these processes. However, urban-scale CTMs often neglect nucleation (by assuming that condensation prevails under conditions with moderate to high PM concentrations) and coagulation (which is typically slow compared to other processes), but regional-scale CTMs with their larger domains and longer transport times may consider these two processes. As described in this section and in Sect. 3.10 and 3.11.1, MACH represents all of the above physical processes as well as several aerosol chemical processes that are described in Sect. 3.5 to 3.7.

3.3.1 Nucleation scheme

Nucleation refers to the creation of new particles by the formation of molecular clusters by low-volatility gases such as sulfuric acid vapour H₂SO₄(g) and some organic species. Such nucleation-mode particles are tiny, with diameters less than 3 nm (Semeniuk and Dastoor, 2018). MACH uses a parameterization of the nucleation process based on classical nucleation theory for sulfuric acid and water binary homogeneous nucleation as proposed by Kulmala et al. (1998) and adopted in CAM (Gong et al., 2003). The expression for particle nucleation rate is dependent on sulfuric acid vapour abundance, temperature, and relative humidity (RH), and the particle mass that is created is assigned to the smallest PM size bin, which for the RAQDPS023 is the PM_2.5 size bin.

3.3.2 Condensation/evaporation scheme

Aerosol particles will grow in size due to the condensation of gas-phase species onto these particles but will shrink in size due to the evaporation of semi-volatile gases from these particles. A modified Fuchs-Sutugin equation is used to calculate the condensation rate for H₂SO₄(g) onto aerosol particles (Fuchs and Sutugin, 1971; Hegg, 1990; see Eq. (A14) of Gong et al., 2003). The expression for the condensation rate to a single particle is linearly dependent on particle diameter and H₂SO₄(g) ambient vapour pressure, where the surface vapour pressure of H₂SO₄ at the particle surface is very small and is assumed to be zero so that H₂SO₄ mass transfer is one-way from the gas phase to the particle phase (Russell et al., 1994). The overall condensation rate to a given particle size bin also depends on the particle number concentration for that bin (N_p,i). This quantity is calculated as the total dry volume of the eight PM chemical components in that bin (i.e., excluding aerosol water) divided by the volume of a spherical particle whose diameter is the arithmetic mean of the size boundaries of that bin. Note that condensation to multiple size bins will occur simultaneously. The condensation process does not change the total particle number density, but it does increase the mass of individual particles. As a consequence, some particles in size bin i may grow in size enough to move to a larger size bin. The treatment of such particles by size “rebinning” is described in Sect. 3.8.

One additional complication is that nucleation and condensation processes compete for H₂SO₄(g). Following the procedure used in CAM (Gong et al., 2003), the division of H₂SO₄(g) between these two processes for each time step is calculated by solving the time rate of change equation of H₂SO₄(g) accounting for both production of H₂SO₄(g) and removal by nucleation and condensation. As noted by Morawska et al. (1999), however, most of the time removal by condensation will dominate removal by nucleation.

The condensation of those organic gases that partition to the particle phase to form SOM (Sect. 3.7) is treated very similarly to the condensation of H₂SO₄(g). The relative condensation rate of H₂SO₄(g) to each PM size bin is used as a proxy for the relative condensation rates of these organic gases to each PM size bin.

3.3.3 Coagulation scheme

Coagulation is the process by which two or more aerosol particles can combine through collision. Coagulation acts to reduce the number concentration of aerosol particles, especially small particles, and hence serves as a control on the evolution of the number concentration of nucleation-mode and directly-emitted accumulation-mode particles. Four distinct physical processes may result in coagulation of particles: Brownian motion; gravitational collection; turbulent inertial motion; and turbulent shear (e.g., Jacobson, 1999; Seinfeld and Pandis, 2016).

MACH follows CAM by using a semi-implicit numerical solution of the general coagulation equation (Jacobson et al., 1994) to compute the coagulation rate and intersectional transfer of aerosol particles (Gong et al., 2003). This scheme conserves bulk particle volume for any time step but results in the transfer of particle volume from smaller to larger size bins. Since MACH assumes that each PM size bin is internally mixed, the mass concentration change of each chemical component can be computed from the calculated volume change of any size bin due to coagulation. Note that the coagulation coefficient is dependent on the “wet”, or real, size of aerosol particles, not the dry size (see Sect. 3.3.5).

3.3.4 Gravitational settling and dry deposition schemes

Particles in the atmosphere undergo gravitational settling (or sedimentation) towards the Earth's surface under the influence of gravity. The downward vertical speed is called the terminal settling velocity and depends on a particle's density and the square of its wet diameter (e.g., Seinfeld and Pandis, 2016). Particles close to the Earth's surface can then be removed from the atmosphere by dry deposition, a complicated process that depends on near-surface meteorological conditions, the nature of the underlying surface, and particle density and size (e.g., Slinn, 1977; Ruijgrok et al., 1995).

Size-resolved particle dry deposition in MACH is represented using a scheme proposed by Zhang et al. (2001), which accounts for a number of relevant processes: turbulent transfer; Brownian diffusion; impaction; interception; gravitational settling; and particle rebound. Calculated particle dry deposition velocities have a strong dependence on particle size with a minimum value at particle diameters of about 1 µm and larger values for both smaller and larger particle sizes (Gong et al., 2003). Surface characteristics for 15 different land-use categories are considered for this calculation (e.g., Zhang et al., 2001; Table S7 of Makar et al., 2018b). Particle dry deposition is only applied at the lowest model level, where the sum of dry deposition velocity and terminal settling velocity is used, whereas gravitational settling occurs throughout the vertical column.

A revised numerical solution for the gravitational settling process was introduced in the RAQDPS023 based on a semi-Lagrangian advection approach. Vertical back trajectories were calculated from the settling and deposition velocities for each particle size, and then mass-conservative interpolation was used to determine the new vertical concentration profile and deposition flux to the surface (Makar et al., 2018a). This new scheme replaced a previous numerical scheme based on an analytical exponential decay rate and better accounts for the model's decreasing vertical grid spacing close to the Earth's surface (Sect. 2.2). The result was a reduction of the rate of removal by downward particle advection. As a consequence, near-surface PM concentrations increased by as much as 0.5 µg m⁻³ compared to those predicted using the old solution scheme (Moran et al., 2021b).

Lastly, to address the strong dependence of particle dry deposition and sedimentation on particle size (e.g., Zhang et al., 2001; Gong et al., 2003), the numerical solution of these processes was modified for the simplified two-bin description of the PM size distribution early in the development of the RAQDPS to reduce numerical errors. The two PM size bins were each subdivided into six sub-bins. Settling velocities were then estimated for each sub-bin based on a third-order polynomial fit to a typical log₁₀(settling velocity) as a function of log₁₀(radius) that was determined using a stand-alone version of the code. Dry deposition velocities were also calculated for each sub-bin based on a weighted average of third-order polynomial fits to a typical log₁₀(dry deposition velocity) as a function of log₁₀(radius) for non-urban land surfaces and urban land surfaces (which has higher values). Revised values of settling velocity and dry deposition velocity for the two size bins were then calculated as logarithms of the sum of the exponentials of the sub-bin values. The revised values of the settling velocities for PM_2.5 and PM_cf based on the sub-bin division were roughly 30 % smaller and 100 % larger, respectively. The revised values of the dry deposition velocities for PM_2.5 and PM_cf, on the other hand, were roughly 100 % larger and 50 % smaller. Analogous modifications to the numerical solution of three other aerosol processes (condensational growth, aqueous-phase chemistry, and inorganic heteorogeneous chemistry) to address intersectional mass transfer in the simplified two-bin configuration are described in Sect. 3.8.

3.3.5 Hygroscopic growth and aerosol activation schemes

The hygroscopic growth of mixed aerosol particles at subsaturation conditions (i.e., RH < 1), which is sometimes referred to as swelling, can result in substantial aerosol water component (WA) values. For example, Tsyro (2005) has suggested that for T=20 °C and RH = 50 % aerosol water may contribute 20 %–35 % of annual mean PM_2.5 concentrations in Europe. Nguyen et al. (2016) have shown based on field study measurements from around the world that aerosol water is ubiquitous and at rural locations contributed on average 3 µg m⁻³ and 35 % of PM₁ total mass. In fact, at high RH values the aerosol water component can dominate aerosol particle composition and mass (e.g., Tang and Munkelwitz, 1994; Nguyen et al., 2016; Widziewicz-Rzońca and Tytła, 2020).

The RAQDPS23 calculates aerosol water using a mixing rule for soluble components (SU, NI, AM, POM, SOM, SS), whereas EC and CM are assumed to be hydrophobic (Gong et al., 2003, Appendix A1). This approach follows that of Hänel (1976). Note that the aerosol water component WA is also diagnosed independently by the HETV scheme (Sect. 3.6).

Under supersaturated conditions (i.e., RH > 1) hygroscopic aerosol particles can grow rapidly by water-vapour condensation to form cloud droplets, whose diameters typically fall in the 10–200 µm range vs. 500–8000 µm for raindrops (e.g., Jacobson, 1999). This process is called aerosol activation but may also be referred to as nucleation scavenging or droplet nucleation, and aerosol particles that can be activated are referred to as cloud condensation nuclei (CCN). Aerosol particle activation is important because it may result in modification of the PM mass and size distributions within and below clouds by aqueous-phase chemistry and wet removal (Sect. 3.5 and 3.10). This trio of processes is referred to collectively as cloud processing. However, whether a CCN particle is activated depends on the level of supersaturation and the particle's size and physical and chemical properties (e.g., Gong et al., 2011).

Two different parameterizations of aerosol particle activation are available in MACH (Gong et al., 2015). For the RAQDPS023, MACH uses the empirical aerosol activation scheme of Jones et al. (1994) to determine the cloud droplet number concentration N_d in a cloud as a function of the particle number concentration N_p. This scheme permits the calculation of a critical particle radius above which all aerosol particles are assumed to be activated (Gong et al., 2006). The portion of aerosol particles incorporated in these cloud droplets is then determined by adding particles from the largest size bins in descending order of size until N_d is reached (Gong et al., 2006, 2015). Next, the bulk cloud liquid water content from GEM (Table 2) is distributed evenly to all activated aerosol particles, after which the cloud droplet size associated with each PM size bin can be determined (Gong et al., 2006). Information about the activated aerosol particles, including their wet diameter, is then available for use by other aerosol process parameterizations, including aqueous-phase chemistry and in-cloud scavenging. Note that one consequence of this scheme is that only some of the particles in a size bin may be activated when the critical particle radius for activation falls between the lower and upper boundaries of a size bin. For the RAQDPS023 with only two size bins, all of the particles in the PM_cf size bin (i.e., N_p,2) must be activated before any particles in the fine size bin can be activated.

3.4 Gas-phase chemistry scheme

Numerous gas-phase chemical reactions occur in the atmosphere (e.g., Jenkin et al., 2015), and some important gas-phase species such as O₃ and hydrogen peroxide (H₂O₂) are reaction products only and are not directly emitted to the atmosphere. It is thus necessary to represent these reactions in a CTM (e.g., Dodge, 2000), but it is also challenging. In particular, volatile organic compounds (VOCs) play a key role in determining the oxidative state of the atmosphere because they may react with three important oxidants: hydroxyl radical (OH); nitrate radical (NO₃); and O₃ (e.g., Atkinson, 1990; Seinfeld and Pandis, 2016). As a consequence VOCs influence photochemistry, PM formation, and acid deposition. However, thousands of VOC species with varying chemical and physical properties, including OH reactivity, vapour pressure, and solubility, are emitted to the atmosphere (e.g., Makar et al., 2003b), and it is not computationally feasible for three-dimensional AQ models to consider the chemical reactions of all of these individual species. Instead, a number of parameterized gas-phase chemistry mechanisms have been developed that consider reduced-form or condensed representations of the full atmospheric chemical system by “lumping” or aggregating groups of VOC species with similar chemical properties together in order to reduce the number of species and chemical reactions that must be considered (e.g., Carter, 1990; Middleton et al., 1990; Kuhn et al., 1998; Dodge, 2000).

The RAQDPS023 used a modified version of the ADOM-2 gas-phase chemistry mechanism to parameterize tropospheric chemistry (Stockwell and Lurmann, 1989; Pudykiewicz et al., 1997; Moran et al., 1998). The ADOM-2 mechanism, which itself was an update of the ADOM-1 condensed mechanism (Lurmann et al., 1986; Fung et al., 1991), was designed to give an accurate representation of hydrocarbon gas-phase reactivity using a modest number of model VOC species. It employs a “lumped molecule” approach (Dodge, 2000) to define 16 lumped VOC species and eight organic radicals in addition to CH₄ and C₂H₆ plus 21 inorganic species. Table 3 lists these 46 ADOM-2 gas-phase chemical species and some of their properties, and Table 4 lists the ADOM-2 mechanism's 114 associated chemical reactions, including 16 photolytic reactions that depend on sunlight.

Table 3 List of ADOM-2 gas-phase model chemical species and their physical and model properties. GMM stands for gram molar mass and DDF stands for dry deposition scaling factor (see Sect. 3.9).

Note from Table 3 that the abundances of 42 of these model gas-phase species were forecast and advected while the abundances of four model species (CH₄, C₂H6, O₂, M) had specified, time-invariant vertical profiles and were not advected or modified by gas-phase chemistry but were considered as reactants in some of the ADOM-2 reactions (Table 4), as was H2O, the water vapour field from GEM (huplus; see Table 2). Eighteen of the advected gas-phase species were organic species, of which 11 were emitted along with eight inorganic species (Table 3). Only two of the ADOM-2 VOC species that were emitted were individual or unlumped species, namely formaldehyde (HCHO) and isoprene (ISOP; C₅H₈). The other nine were lumped species, and this high degree of lumping imposes a limitation for evaluating model predictions of VOC species (e.g., Stroud et al., 2008; Moran et al., 2026).

Table 4 List of ADOM-2 gas-phase mechanism chemical reactions and reaction rate coefficients in KPP-compliant input format (Damian et al., 2002; Sandu and Sander, 2006). Chemical reactions are listed in the left-hand column, and include first-, second-, and third-order reactions. All of the chemical species are defined in Table 3, except for hυ, which indicates a photon of light of frequency υ, and DUMMY, which indicates a non-reacting species. The reaction rate coefficients listed in the right-hand column have units of s⁻¹, cm³ mol⁻¹ s⁻¹, and cm⁶ mol⁻² s⁻¹ for first-, second-, and third-order reactions, respectively. SUN is normalized sun intensity and ranges from 0 at night and 1 at noon (but is replaced during simulations by model-calculated J values); TEMP corresponds to ambient atmospheric dry-bulb temperature (K); the ARR2 function describes the reaction rate coefficient in Arrhenius form, where the first argument is the preexponential factor and the second argument provides the $E/R$ value; the TYPE5 function calculates third-order, pressure- and temperature-dependent reaction rate coefficients based on Troe theory (e.g., Seinfeld and Pandis, 2016), where the first argument specifies the F factor, the second and third arguments specify the high-pressure limiting rate coefficient, and the fourth and fifth arguments specify the low-pressure limiting rate coefficient; RCONST(n) is the reaction rate coefficient for the nth reaction; and CFACTOR is a units conversion factor.

Note that Reaction (R21) represents a two-step heterogeneous reaction on a water surface and that Reactions (R36) and (R42) have been combined with Reactions (R35) and (R41), respectively.

The calculation of lumped ADOM-2 VOC emissions begins with the assignment of individual emitted VOC species to one of the 32 VOC categories proposed by Middleton et al. (1990) to aggregate the 1980 National Acid Precipitation Assessment Program (NAPAP) anthropogenic VOC emissions inventory or to three additional biogenic categories (Isoprene, Alpha- and Beta-Pinene, and Other Monoterpenes) suggested by Makar et al. (2003b). These 35 VOC categories are then further lumped to the 11 emitted ADOM-2 VOC species using mole-based reactivity weights (see Table 5). Most of the ADOM-2 VOC species represent tens, hundreds, or even thousands of individual VOC species. For example, the ADOM-2 C3H8 lumped VOC species includes propane (C₃H₈) as would be expected, but it also includes benzene, acetylene, and some other lower-reactivity species (e.g., alcohols, ethers, and esters) that are assigned to NAPAP category 27 (e.g., Plummer et al., 2001; Makar et al., 2003b). The ETHE lumped VOC species includes ethene (C₃H₈), but it also includes isoprene oxidation products such as methacrolein and methyl vinyl ketone in order to implement a condensed isoprene chemistry mechanism used by the Carbon Bond-IV chemistry mechanism (Gery et al., 1988, 1989; Stockwell and Lurmann, 1989; Stroud et al., 2008). Note that emissions of acetone, organic acids, haloalkenes, and three other NAPAP VOC categories (26, 31, 32) are not included in emissions for the ADOM-2 mechanism due to their low or unknown reactivities (referred to collectively in Table 9 as EOTH). While the omission of these categories is defensible from a reactivity perspective, the neglect of organic acids in particular may remove a potential source of secondary organic aerosol (e.g., Makar et al., 2003b; Fisseha et al., 2004).

Table 5 List of NAPAP VOC classes and properties and mapping to ADOM-2 lumped VOC species. The “Others” classes include alcohols, ethers, alcohol ethers, esters, etc., while OTH includes all low- and non-reactive species.

Rate constants for the inorganic reactions considered by the ADOM-2 mechanism were based largely on the recommendations of DeMore et al. (1987) while those for the organic reactions were based largely on Atkinson et al. (1992). The rates of the 16 photolytic reactions in Table 4 are expressed as the abundance of photoactive species multiplied by the photodisassociation rate coefficient (i.e., J value). The ADOM-2 photodisassociation rate coefficients were calculated from lookup tables of clear-sky actinic flux from Peterson (1976), which were based on the radiative transfer model of Dave (1972) and absorption cross-sections and photodisassociation quantum yields from DeMore et al. (1987) (Stockwell and Lurmann, 1989). The photodisassociation rate coefficients for $J_{{\mathrm{NO}}_{\mathrm{2}}}$ and $J_{{\mathrm{O}}_{\mathrm{3}}\to \mathrm{O}{(}^{\mathrm{1}}\mathrm{D})}$ depend on both altitude and solar zenith angle whereas the other photodisassociation rate coefficients depend only on solar zenith angle and species-dependent scale factors (Kelly et al., 2012). Lastly, to account for the effect of clouds on clear-sky photolysis rate constants, the total cloud fraction field f_cld (ftot in Table 2), which accounts for both resolved and SGS clouds in a model layer, and cloud liquid water content field (lwc in Table 2) predicted by GEM were used to scale the precalculated clear-sky J values following an algorithm from Chang et al. (1987). More details are given in Majdzadeh et al. (2022). Note that the performance of the ADOM-2 gas-phase mechanism was compared to a number of other gas-phase chemistry mechanisms by Kuhn et al. (1998), and its predictions were found to be close to the median for the nine mechanisms tested.

To integrate the chemistry tracer conservation-of-mass equations in time, the Young and Boris (1977) predictor-corrector method was used to solve the gas-phase chemistry step of the process operator splitting sequence (see Sect. 4.1). This asymptotic algorithm, which typically must use much smaller integration time steps than the overall MACH chemistry time step of 900 s due to the stiffness of the ADOM-2 chemical system, is integrated over the MACH chemistry time step.

The treatment of stratospheric chemistry in the RAQDPS023 was much simpler than that for tropospheric chemistry. The LINOZ scheme for ozone stratospheric chemistry (McLinden et al., 2000) was used to forecast ozone concentrations while concentrations of other species were only affected by dynamics and physics. The definition of the stratosphere used to apply the LINOZ scheme in each grid column was either those vertical levels located above 100 hPa or with specific humidity less than 10 ppmv. The ADOM-2 scheme was used below these levels. The LINOZ scheme was also implemented in the GDPS 8.0.0 to support the data assimilation of satellite-based ozone measurements (CMC-GDPS-8.0.0, 2021a).

3.5 Aqueous-phase chemistry scheme

Cloud chemistry, that is, aqueous-phase chemistry in cloud water, is an important pathway for the conversion of SO₂ to SO${}_{\mathrm{4}}^{=}$ (e.g., Barth et al., 2000; Gong et al., 2011). As well as referring to chemical reactions amongst various species in dilute aqueous solution in the cloud droplets, the term aqueous-phase chemistry usually encompasses two related processes: mass transfer of species between the gas phase and aqueous phase (absorption/condensation) inside clouds and the dissociation/ionisation in cloud water of certain dissolved species. Note that consideration of these processes is necessarily restricted to activated aerosol particles (Sect. 3.3.5).

MACH employs a slightly adapted version of the ADOM aqueous-phase chemistry mechanism, which considers 12 gas-phase and 13 aqueous-phase species together with 25 reactions, to represent these processes (Young and Lurmann, 1984; Fung et al., 1991). Note from Table 4 that one ADOM-2 mechanism species, NH₃, does not participate in any gas-phase reactions, but NH₃ gas participates in aqueous-phase and inorganic heterogeneous reactions (see also Sect. 3.6). Table 6 lists the 25 gas and aqueous species and Table 7 lists the 25 reactions that are considered by this condensed mechanism. Reactions (R19)–(R22) describe the four pathways for SO₂ oxidation (i.e., S(IV) → S(VI)) included in the mechanism: these reactions involve three aqueous-phase oxidants (H₂O₂, O₃, ROOH) and oxygen catalyzed by iron and manganese (e.g., Seinfeld and Pandis, 2016). Fourteen of the mechanism reactions in Table 7 describe the mass transfer of soluble gases as a reversible diffusion process constrained by chemical equilibrium (i.e., forward-backward reaction pairs R2–R3, R4–R5, R6–R7, R8–R9, R10–R11, R12–R13, and R15–R16 for gas-phase SO₂, O₃, H₂O₂, HNO₃, ROOH, NH₃, and CO₂, respectively). Moreover, reaction pairs (R2)–(R3), (R8)–(R9), and (R12)–(R13) describe both mass transfer from the gas- to the aqueous phase and subsequent dissociation in cloud water. Aerosol particle scavenging, which is represented by Reactions (R1) and (R23)–(R25), is an irreversible process that can occur by nucleation, Brownian motion, phoretic attachment, or inertial impaction (e.g., Wang et al., 2010), but only nucleation (i.e., aerosol activation) is currently considered for in-cloud scavenging. Inorganic particle components SU, NI, and AM enter cloud droplets through this process. The diffusion coefficients needed to describe the mass transfer process are determined using the Fuchs and Sutugin (1971) formulation (Young and Lurmann, 1984).

Table 6 List of ADOM aqueous-phase mechanism chemical species and phases (G = gas, P = particle, A = aqueous).

Table 7 List of ADOM aqueous-phase mechanism chemical reactions, where gas-phase abundance units are ppmv and aqueous-phase concentration units are mol L⁻¹ (see Young and Lurmann, 1984).

L_w= volumetric liquid water fraction in air (m³ H₂O per m³ air). R= universal gas constant (8.205×10⁴ ppm L mol⁻¹ K⁻¹). T= dry-bulb temperature (K). k= equilibrium constant for aqueous-phase reactions. B1 =L_wRT (conversion factor: gas-phase to aqueous-phase units). B2 = (L_wRT)⁻¹ (conversion factor: aqueous-phase to gas-phase units). B3, B4, and B5 are defined by the emissions data.

To reduce computer time, following Gong et al. (2006) this chemistry mechanism is implemented in MACH as a “bulk” (i.e., non-size-resolved) process that occurs in a generic cloud droplet whose size is determined from the bulk or total cloud liquid water content (LWC) supplied from GEM (Table 2) and cloud droplet number concentration N_d (Sect. 3.3.5). In addition, the GEM total cloud fraction field f_cld is used to convert grid-average values to “in-cloud” values before the cloud chemistry operator is applied.

Like the gas-phase mechanism the aqueous-phase mechanism is solved for a MACH chemistry time step using the computationally efficient Young and Boris (1977) predictor–corrector algorithm with a number of smaller integration time steps. At the end of the aqueous-phase chemistry integration step, the gridded gas-phase concentrations of SO₂, O₃, H₂O₂, HNO₃, ROOH, and NH₃ are updated to account for their in-cloud depletion and the resulting bulk mass increments of three dissolved inorganic aerosol components (SU, NI, AM) are redistributed across activated PM size bins in air, effectively returning particle concentrations to clear air conditions. This redistribution is done by using the ratios of LWC in each activated (or partially activated) size bin to the total LWC, which implies that aqueous-phase chemistry is a volume-controlled process (Gong et al., 2006). The LWC in each bin depends on the number of activated aerosol particles in each bin and is obtained by distributing bulk LWC to all activated aerosol particles evenly (Gong et al., 2006, 2011). This mass redistribution step at the end of the aqueous-phase chemistry operator is then followed by a “rebinning” step that is needed to maintain the fixed size bin structure (Sect. 3.8). Note, however, that the resulting ambient concentrations of inorganic particle components NI and AM may not be in equilibrium with the gas phase so the inorganic heterogeneous chemistry operator is called next.

3.6 Inorganic heterogeneous chemistry scheme

Ambient sulfuric and nitric acid vapours and ammonia gas can partition to atmospheric particles to form one or more inorganic salts, and the resulting SU, NI, and AM aerosol components can make up a significant fraction of PM_2.5 mass (e.g., Brook and Dann, 1999; Malm et al., 2004). Depending on the relative fractions of these three gas-phase species, four different salts may be formed (e.g., Ansari and Pandis, 1998): ammonium sulfate ((NH₄)₂SO₄), ammonium bisulfate (NH₄HSO₄), letovicite ((NH₄)₃H(SO₄)₂), or ammonium nitrate (NH₄NO₃). Because H₂SO₄(g) has a very low vapour pressure (Sect. 3.3.2), it is a reasonable approximation to assume that it resides completely in the aerosol phase (Nenes et al., 1998), whereas nitric acid vapour (HNO₃(g)) and ammonia gas (NH₃(g)) exist in a reversible equilibrium between the gas and particle phases. Knowledge of inorganic heterogeneous chemistry thus allows the secondary inorganic component of PM to be determined.

Inorganic heterogeneous chemistry was parameterized in the RAQDPS023 version of MACH with the computationally efficient HETV (HETerogeneous Vectorized) scheme (Makar et al., 2003a), which is based on a subset of the ISORROPIA thermodynamic equilibrium algorithms (Nenes et al., 1998, 1999). In ISORROPIA the full system of inorganic equilibrium equations is broken up into subsystems, each corresponding to fixed ranges of RH and the ratio of total ammonia to total sulfate, where total ammonia includes both gas and particle forms in molar units. These ranges constrain which inorganic species (gases, ions and salts) may be present in aerosol particles. HETV considers 12 such range-based subsystems or “cases”. The multicomponent activity coefficient calculations that are required for non-ideal, high-concentration solutions use the formula of Kusik and Meissner (1978) as described by Kim et al. (1993) but are approximated using higher-order Taylor series to reduce calculation time (Makar et al., 2003a). The metastable state assumption, namely that some particle liquid water is always present even at low RH, was also made to reduce the number of cases that must be considered (e.g., Rood et al., 1989; Miller et al., 2024).

Similar to the implementation of the aqueous-phase chemistry parameterization (Sect. 3.5), HETV assumes a bulk equilibrium in the particulate phase over all particle sizes to save computer time. At the end of the inorganic-heterogeneous-chemistry integration step, gridded concentrations of HNO₃(g) and NH₃(g) are updated to account for inorganic gas-particle partitioning during the time step, where the updated values may have either increased or decreased. If the bulk mass increments of the three inorganic aerosol components (SU, NI, AM) are positive, then these increments are redistributed across the PM size bins. This redistribution is carried out by using the ratios of gas-phase condensation rates to all particles in a size bin based on the modified Fuchs-Sutugin equation (Hegg, 1990; Makar et al., 1998; Gong et al., 2003) and the total condensation rate summed over all size bins. If the bulk mass increments are negative, though, then mass is removed from each size bin starting with the largest based on the same condensation rate ratios. Note that it is possible that for some bins the mass that is available to be removed may be less than the pro-rated increment, in which case a second removal iteration may be required. This redistribution step is then followed by a mass-conserving “rebinning” step that is needed to maintain the fixed size bin structure (Sect. 3.8).

3.7 Secondary organic aerosol formation scheme

The term “secondary organic aerosol” (SOA; equivalent to SOM) refers to organic PM that results from the oxidation of gas-phase VOC precursor species followed by condensation to the particle phase (e.g., Hallquist et al., 2009; Jimenez et al., 2009). The RAQDPS023 used the Instantaneous secondary organic Aerosol Yield (IAY) scheme described by Jiang (2003, 2004, 2005) to represent SOA formation from five model VOC species. Four of these are lumped VOC species: ALKA (long-chain anthropogenic and biogenic alkanes); ALKE (long-chain anthropogenic and biogenic alkenes); AROM (multi-substituted aromatics); and TOLU (toluene and mono-substituted aromatics) (Stroud et al., 2008). The fifth model VOC species is ISOP (isoprene), which is not lumped and is dominated by biogenic emission sources (although there are also minor anthropogenic sources). The total SOA contribution from ALKE oxidation is further divided into contributions from anthropogenic ALKE and from biogenic α-pinene and β-pinene (derived from BEIS monoterpene emissions; see Sect. 3.11.3).

The Jiang scheme is based on a two-condensable-product fit to chamber data (Pankow, 1994; Griffin et al., 1999), where initial products are assumed to be converted rapidly to non-volatile organic PM. Values of the IAY parameters that were used by the RAQDPS023 are listed in Table 8, where α_1i and α_2i are the mass-based stoichiometric coefficients (unitless) and K_om,1i and K_om,2i are the partitioning coefficients (m⁻³ µg) of the ith condensable species CS_i, respectively (Odum et al., 1996; Jiang, 2003). The α_i and K_om,i values for ALKA and anthropogenic ALKE are taken from Griffin et al. (1999), those for AROM, α-pinene, and β-pinene are taken from Jiang (2003), the TOLU values come from Odum et al. (1997), and the ISOP values are taken from Barsanti et al. (2013).

The amount of SOA product available for condensation at each MACH chemistry time step is the sum of IAY_i⋅ CS_i over the seven VOC condensable species (CS) values listed in Table 8, where IAY_i is a dimensionless fraction and CS_i is the reduction due to oxidation of ALKA, ALKE, AROM, TOLU, and ISOP MMRs over a MACH time step as calculated by the gas-phase chemistry operator (Sect. 3.4). The condensation rate of SOA products to aerosol particles is described by a modified Fuchs-Sutugin equation (see Eq. (A14) of Gong et al., 2003), similar to the treatment for H₂SO₄ condensation (Sect. 3.3.2), and the relative condensation rates of H₂SO₄(g) to each PM size bin are used as a proxy for the relative condensation rates of these SOA organic gases. Lastly, like the H₂SO₄ condensation step the SOA condensation step must be followed by a “rebinning” step to maintain the fixed PM size bin structure as described next.

Table 8 List of IAY parameter values for seven, SOA-forming, ADOM-2 mechanism VOC species. α_i is the mass-based stoichiometric coefficient (unitless) and K_om,i is the partitioning coefficient (m⁻³ µg) of the ith condensable species.

3.8 Intersectional mass transfer schemes

As already noted, one complication of using a fixed sectional representation of the PM size distribution is that particles belonging to a particular size bin may grow too large for that bin due to the addition of mass from condensation or from aqueous-phase or heterogeneous chemistry (e.g., Jacobson, 1999; Gong et al., 2006). This problem can be addressed by applying a rebinning operator that partitions the enlarged particles between two adjacent size bins in a mass- and number-conserving manner in order to maintain the fixed bin structure. Jacobson (1999) refers to this rebinning approach as a quasi-stationary size structure for the PM size distribution. Interestingly, fixed modal representations of the PM size distribution suffer from an analogous issue that has an analogous solution, referred to as mode merging (e.g., Binkowski and Shankar, 1995; Binkowski and Roselle, 2003).

MACH implements a mass-conservative rebinning procedure each time step after each of three processes: condensational growth (Sect. 3.3.2 and 3.7); aqueous-phase chemistry (Sect. 3.5); and inorganic heterogeneous chemistry (Sect. 3.6). A slightly different approach to rebinning is used for each of these three processes. For condensational growth the rebinning calculation is treated as mathematically equivalent to two-particle coagulation (Gong et al., 2003). The volume of a particle in the size bin that underwent condensation (the “donor” bin) is considered to be one colliding particle while the volume of the condensable species per particle is considered to be the other colliding particle. The same formulation used by Jacobson et al. (1994, Eq. 13) to partition the volume of an intermediate particle of two coagulating particles into two model size bins (the “receiver” bins) is then used (Gong et al., 2003). The two volume partitioning factors that are calculated depend on the particle volumes of the donor and receiver size bins before condensation occurred and the particle volume for the donor bin after condensation occurred. The sum of the two partitioning factors is unity (to conserve the allocated volume), and the smaller of the two adjacent receiver size bins may be the same as the donor size bin but may also be a larger size bin (if n>2). The choice of the receiver size bins depends on the mean particle diameter diagnosed for the new total mass for the donor size bin mass after condensation using the initial donor-bin particle number. The rebinning calculations are applied to each size bin beginning with the smallest one. Note that all PM chemical components are affected by rebinning due to the assumption of internal mixing even though condensational particle growth is due only to one or two chemical components (i.e., sulfuric acid vapour or condensable organic gases).

For activated particles that have grown due to the production of particle sulfate by aqueous-phase chemistry, rebinning is performed after the redistribution step at the end of the aqueous-phase chemistry operator. A similar equivalence to two-particle coagulation is invoked, with each activated particle size bin corresponding to a donor bin and the size-bin mass gain due to cloud chemistry corresponding to the second particle, but the partitioning-factor pairs for the receiver bins are calculated with the mass-based method described by Gong et al. (2006; Eq. 3), which also conserves particle number. Rebinning is restricted to size bins with activated particles, where the rebinning calculations are applied beginning with the smallest activated or partially activated bin and then moving to successively larger activated size bins but are based on non-activated particle properties.

The inorganic heterogeneous operator is called immediately after the aqueous-phase chemistry operator (see Sect. 4.1) because the updated PM composition resulting from the latter may not be in stable equilibrium with the gas phase However, once HETV determines a new gas-particle partitioning equilibrium, some particles may have gained mass or lost mass and hence may need to be assigned to a different size bin. The mass-based method of Gong et al. (2006) is also used for rebinning here but with two adjustments. First, activation is not relevant so in the case of mass gain rebinning is applied to all particle size bins beginning with the smallest size bin and moving to successively larger bins. Second, since particles may have lost mass as well as gained mass, the two-particle coagulation analogy is generalized to account for a “donor” particle with a negative mass gain. For such negative mass increments the rebinning moves mass to a pair of smaller receiver size bins beginning with the largest size bin and moving to successively smaller bins.

The above descriptions are general and also apply to the 12-bin version of GEM-MACH. For the RAQDPS023, which only considers two PM size bins, the two adjacent receiver size bins can only be the PM_2.5 and PM_cf size bins and the rebinning step can only transfer particle mass from the PM_2.5 size bin to the PM_cf size bin. During development of the original GEM-MACH-based RAQDPS, it was found that the rebinning operator was transferring too much mass from the PM_2.5 size bin to the PM_cf size bin. This problem was not apparent for the 12-bin version, suggesting that the very coarse representation of the PM size distribution was introducing numerical errors. Conceptually if we consider that it is those particles whose diameters are already close to the upper end of a size bin that are most likely to grow into the next larger size bin after mass is added by condensation or other processes, then an assumption that particles of any size from the ultrafine to 2.5 µm might increase enough in size in one time step to join the PM_cf size bin is unrealistic.

The solution adopted for the two-bin version to avoid this assumption was to subdivide each size bin into a number of sub-bins, allocate the mass increase for each size bin across the sub-bins, apply the rebinning operator to the individual sub-bins, and then aggregate the sub-bin particles back to the original size bins. In addition to the choice of the number of sub-bins, however, the subdivision and allocation steps could be performed in a number of ways, such as a linear vs. a logarithmic subdivision of size bins and allocation of the mass increment to a uniform distribution vs. an assumed log-normal distribution.

After some experimentation a linear subdivision of each size bin into six equal-sized sub-bins by radius was adopted with the mass of each chemical component for the size bin divided equally between the six sub-bins. For condensational growth the allocation of size-bin gained mass for a time step across the sub-bins was based on the relative fraction of total particle area by sub-bin for the size bin. For cloud chemistry the allocation step was done in three steps. First, aerosol number concentration was allocated to sub-bins based on the equal distribution of particle mass to the sub-bins. Next, cloud liquid water was allocated to each sub-bin based on the relative fraction of total particle number by sub-bin for the size bin. And third, the relative fraction of cloud liquid water by sub-bin for the size bin was used to allocate the mass increment to each sub-bin. Finally, for inorganic heterogeneous chemistry the allocation step for aerosol particle mass changes followed the approach for condensational growth: the size-bin mass increment for a time step was divided equally between the six sub-bins, but with the additional possibility that the size-bin mass increment could be negative, that is, a mass decrement, due to volatilization of the nitrate and ammonium aerosol components. The net result of introducing these modifications was to reduce the intersectional transfer of mass after these process operators from the PM_2.5 bin to the PM_cf bin.

3.9 Gas-phase dry deposition scheme

Many atmospheric gases, like aerosol particles (Sect. 3.3.4), are removed at the Earth's surface by dry deposition after downward transport by air motions (though gravitational settling plays no role, unlike the particle case). The rate of removal by dry deposition is determined in part by the chemical and physical properties of each individual gas, including its reactivity, solubility, and molecular diffusivity. The rate of removal is also affected by surface properties, including vegetation characteristics in the case of removal through plant stomata and by other vegetative surfaces, including wetted surfaces (e.g., Wesely and Hicks, 2000; Zhang et al., 2002b; Clifton et al., 2020, 2023; Galmarini et al., 2021).

To represent the dry deposition of gas-phase species MACH uses a variant of Wesely's resistance-analogy parameterization (Wesely, 1989), which is based on the product of the concentration of a gas-phase species at the Earth's surface and its dry deposition velocity. Following an analogy to electrical circuits, the dry-deposition velocity of a species is given by the inverse of total “resistance”, which consists of the sum of three components: aerodynamic resistance, quasi-laminar sublayer resistance, and surface resistance. Aerodynamic resistance R_a is a function of only near-surface micrometeorological conditions and the roughness characteristics of the underlying surface, so it is independent of the chemical species. Quasi-laminar sublayer resistance R_b is a function of friction velocity and the molecular diffusivity of each chemical species. The surface resistance R_c is the most important and complex resistance for most gas-phase chemical species, and it effectively determines whether a particular species is dry deposited or not. Soil resistance, soil wetness, and snow cover must be considered for all land surfaces in calculating R_c, with additional resistances (stomatal, mesophyll, cuticle, in-canopy, exposed surfaces) required in the case of vegetated surfaces (Makar et al., 2018b).

To account for the different dry deposition rates of gas-phase species, MACH considers two master species, SO₂ and O₃, which behave quite differently for wet surfaces (Zhang et al., 2002a, b). Weighted averages of the dry deposition velocities for these two master species are then used to scale dry deposition velocities for other dry depositing species (Zhang et al., 2002a; Makar et al., 2018b). Table 3 lists the 17 ADOM-2 gas-phase species that are dry deposited and the master-species dry deposition scaling factors (DDF) that are assumed for each. Stomatal resistance is parameterized using a one-big-leaf approach, where the effects of incoming solar radiation, temperature, and humidity deficit are empirically represented but the increase of shading due to increasing leaf area index (LAI) is neglected (Makar et al., 2018b).

One limitation of the gas-phase dry deposition scheme, however, is its simplistic treatment of vegetation phenology and seasons in both time and space: only five phenologically-based seasons are considered, which vary by month and 5° latitude bands without any dependence on longitude or elevation (Brook et al., 1999a). These seasons are used in lookup tables to specify seasonal values of a number of required parameters, including soil, cuticle, canopy, and minimum stomatal resistances, for 15 land-use types (Makar et al., 2018b). For some months, however, this treatment can result in zonal bands being visible in surface concentration fields due to seasonal differences between latitude bands. Two removal processes for SO₂ were also recently identified as missing from the RAQDPS023: the soil-wetness and cuticle-wetness gas-phase dry deposition pathways (e.g., Wesely, 1989; Galmarini et al., 2021; Clifton et al., 2023). Further details about the scheme used by MACH to model the dry deposition of gases may be found in Makar et al. (2018b), and evaluation results for O₃ deposition for this scheme and for comparable schemes used by other CTMs can be found in three recent papers (Clifton et al., 2023; Hogrefe et al., 2025; Kioutsioukis et al., 2025).

3.10 Wet deposition schemes

Wet deposition, which is also referred to as wet removal or precipitation scavenging, consists of multiple processes related to tracer scavenging and transport by hydrometeors, which may have the form of raindrops, snowflakes, ice pellets, graupel, or hail (e.g., Slinn, 1984; Seinfeld and Pandis, 2016). Some wet deposition processes act within clouds while other wet deposition processes act below clouds. The former are referred to as in-cloud scavenging or “rain-out” while the latter are referred to as below-cloud scavenging or “wash-out”. “Snow-out” must also be considered, and, unlike raindrops, solid hydrometeors can have many shapes or habits, including plates, needles, dendrites, and columns (e.g., Slinn, 1984; Zhang et al., 2013). As well, the processes responsible for the removal of soluble gases vs. aerosol particles are not identical. But regardless of these complexities, a number of modelling studies have suggested that wet deposition can be the dominant removal pathway for aerosol particles, though not for soluble gas-phase species (e.g., Barrie et al., 2001; Iversen and Seland, 2002; Textor et al., 2006).

As described in Sect. 3.3.5 aerosol particles that are activated in the supersaturated conditions within clouds act as CCN and become cloud droplets. Soluble gases inside clouds may also enter cloud droplets through gas-particle partitioning (Sect. 3.5). The majority of clouds are non-precipitating, but if cloud droplets do form raindrops that then fall from the cloud, both dissolved particle mass and dissolved gases from the cloud layer are also carried downwards by the raindrops. This “rain-out” process is parameterized in MACH using both bulk autoconversion (or precipitation production or cloud-to-rain conversion) rate f_ctr obtained from GEM (ppro in Table 2; see Eq. (5) of Gong et al., 2006) and SGS convective precipitation production from the middle and deep convection implementations of the Kain-Fritsch scheme (Sect. S1.7). The GEM total cloud fraction field f_cld (ftot in Table 2) is then used to convert grid-averaged parameters (e.g., cloud liquid water content, precipitation rates/fluxes) for each model layer to “in-cloud” values (Gong et al., 2006). The contribution of mid-level SGS convective precipitation to precipitation scavenging was added in MACH for the RAQDPS023 (Moran et al., 2021b). Note that the “snow-out” process within cold clouds consisting of ice crystals is not considered explicitly by MACH (Gong et al., 2006).

Falling hydrometeors can also scavenge aerosol particles and soluble gases located below cloud base (e.g., Slinn, 1984; Gong et al., 2006; Wang et al., 2010; Zhang et al., 2013). To calculate the below-cloud scavenging of aerosol particles by liquid or solid precipitation, MACH follows CAM (Gong et al., 1997a, 2003; Gong et al., 2011; Ghahreman et al., 2024) and employs a simplified treatment of the particle-size-dependent precipitation scavenging rate coefficient Λ(d) proposed by Slinn (1977, 1984), where d is the particle diameter. This includes a semi-empirical expression for the collection efficiency E(d,D), where D is a characteristic size of the hydrometeor, that accounts for the three dominant particle collection mechanisms (Brownian diffusion, interception, inertial impaction) and assumes that the retention efficiency is unity (Slinn, 1977, Eq. 10). In the case of rain, D corresponds to the diameter of a spherical raindrop, the particle settling velocity is neglected relative to the raindrop terminal fall speed V_rt, and the precipitation number distribution is assumed to be monodisperse and characterized by a mean raindrop diameter D_m that can be calculated from the precipitation rate p (Slinn, 1977, Eq. 14). V_rt is calculated using a formula from Beard (1976). D_m and V_rt are then used with Eq. (10) of Slinn (1977) to calculate the mean collection efficiency $\overline{E}(d,D_{\mathrm{m}}$), which is used in turn to calculate Λ_r(d), the particle-size-dependent rate coefficient for rain scavenging. Λ_r(d) depends directly on p and $\overline{E}(d,D_{\mathrm{m}}$) and inversely on D_m (Slinn, 1977, Eq. 13).

A similar formulation is used by MACH to parameterize below-cloud scavenging of particles by snow, but the precipitation rate p corresponds to a rainwater equivalent (snow density is assumed to be a factor of 10 less than rain) and the mean collection efficiency $\overline{E}(d,l)$ depends on l, a characteristic dimension of the solid hydrometeors for particle collection that may be different from D, the size of the hydrometeors. The resulting snow scavenging rate coefficient Λ_s(d) depends directly on p and $\overline{E}(d,l)$ and inversely on D_m (Slinn, 1984, Eq. 11.36). The frozen hydrometeor terminal fall speed V_st needed to calculate $\overline{E}(d,l)$ is obtained using the parameterization of Slinn (1984) (see also Gong et al., 1997a). Slinn (1984) suggested D_m values for five different snow habits, of which MACH considers snow scavenging for three different temperature ranges by needle snow (−8 to 0 °C), stellar snow (−25 to −8 °C), and graupel ($<-\mathrm{25}$ °C) (Gong et al., 2011).

Below-cloud scavenging is assumed to be due to a single precipitation phase, either rain or snow. A threshold of the ambient dry-bulb temperature T is used to determine the phase: rain for T>0 °C and snow for T < 0 °C (Ghahreman et al., 2024). To account for wet removal of particles in the partial vertical column below cloud base, a two-dimensional cloud cover field is first calculated from the three-dimensional total cloud fraction field f_cld using a random overlapping algorithm (Mailhot et al., 1998). This two-dimensional cloud cover field is then used to convert grid-averaged parameters (e.g., precipitation rate) into separate below-cloud and clear-sky values and then to recalculate grid-averaged values (Gong et al., 2006).

Six soluble gases (SO₂, HNO₃, NH₃, H₂O₂, ROOH, CO₂) are also removed in MACH by below-cloud precipitation scavenging. The removal of HNO₃, NH₃, H₂O₂, and ROOH is treated as an irreversible process similar to particle scavenging. This is a reasonable assumption when the ambient concentration of the gas-phase species is much greater than its equilibrium concentration at the surface of the raindrop or if the effective Henry's law coefficient for the species is very large (e.g., Seinfeld and Pandis, 2016). Λ_i, the scavenging rate coefficient for rain for the ith gas species, is calculated using Eqs. (6) and (7) of Gong et al. (2006), where the four key parameters are p, D_m, V_rt, and D_i, the gas-phase diffusivity of the ith species. Wet scavenging of SO₂ and CO₂, on the other hand, is treated as a reversible process since their gas-phase and gas-droplet equilibrium concentrations may be closer in magnitude or their effective Henry's law coefficient is not large. The two equilibrium equations given by Eq. (8) of Gong et al. (2006) are used for these two species, where the key parameter is H^′, the effective equilibrium constant that combines absorption and dissociation equilibria. Lastly, scavenging of soluble gases by snow and ice is only considered for HNO₃ and NH₃. This process is assumed to be irreversible and is based on scaling the H₂SO₄ scavenging rate (Karamchandani et al., 1985; Gong et al., 2006).

While wet deposition results from hydrometeors transporting pollutants downwards to the Earth's surface where they are removed, some hydrometeors may evaporate before they reach the ground, in which case the pollutants are released back to the ambient air, though at a lower elevation. This vertical transport process is explicitly included in MACH (see Gong et al., 2006) by using the precipitation evaporation rate predicted by GEM (pevp in Table 2). Pollutants are assumed to be released from evaporating hydrometeors at the same rate as water vapour.

3.11 Chemical emissions inputs

Files of gridded chemical emission fluxes are a key input to CTMs, and they have a direct and large impact on model performance. In effect these emission fluxes constitute a chemical forcing term that transforms the chemical-species mass conservation equations of a CTM from an initial-value problem into an initial-boundary value problem. However, the preparation of these emission fluxes into a model-ready form is a challenging task. For one thing, all sources of emission fluxes must be considered, including both anthropogenic and natural emissions. In addition, extensive pre-processing is required to prepare files of anthropogenic emission fluxes before a CTM simulation can be performed whereas natural emissions are typically not known in advance and must be modelled and predicted during a CTM simulation.

A new set of gridded, hourly anthropogenic emissions flux input files for each month of the year for a typical seven-day week was generated for use by the RAQDPS023 and RAQDPS-FW023 (Moran et al., 2021b; see also Table A1). This emissions data set is referred to as SET4.0.0. This section describes the data sources and the preparation of the SET4.0.0 data set along with the modelling of biomass burning emissions, biogenic emissions, and sea-salt emissions.

3.11.1 Anthropogenic emissions

Information about year-specific anthropogenic emissions for North America for use by CTMs is usually obtained from available annual national emissions inventories (e.g., Miller et al., 2006; Moran et al., 2013; Day et al., 2019; Foley et al., 2023). Such inventories are most often compiled by government agencies. The following three annual national anthropogenic emission inventories were used as the basis for generating the SET4.0.0 emissions input files:

• Projected 2020 Canadian Air Pollutant Emissions Inventory (APEI) based on the 2015 APEI (ECCC, 2016; Sassi et al., 2021);

• Projected 2023 U.S. National Emissions Inventory (NEI) based on version 1 of the U.S. Environmental Protection Agency (EPA) 2016 Emissions Modeling Platform for policy development applications (U.S. EPA, 2019); and

• Projected 2023 Mexican NEI based on 2008 Mexican NEI (U.S. EPA, 2019, p. 168).

Each of these inventories contains annual or monthly emissions of seven criteria air contaminants (CACs), namely SO₂, NO_x, CO, VOC, NH₃, PM_2.5, and PM₁₀, reported for different emission source sectors (e.g., heavy-duty diesel vehicles, residential wood combustion, solvent use, agricultural tillage, etc.) by subnational jurisdiction. Emission source sectors corresponding to 417 Source Classification Codes (SCCs) were contained in the 2020 Canadian inventory vs. 5473 SCCs for the 2023 U.S. inventory and 791 SCCs for the 2023 Mexican inventory, while Canadian emissions were reported for 13 provinces and territories vs. 3249 counties for the U.S. and 2457 counties for Mexico (CMC-RAQDPS-023, 2021). It should be noted that the level of detail or granularity of emissions reporting in each inventory affects the level of detail that can be employed for subsequent emissions processing (see Sect. S2). For example, the county-level reporting of emissions for the U.S. and Mexico vs. province-level reporting for Canada meant in effect that some spatial allocation had already been performed for the U.S. and Mexico. Note too that not all U.S. and Mexican counties were located within the RAQDPS023 model domain, which meant that these emissions were not included in the SET4.0.0 files.

Most emissions inventories are retrospective (i.e., historical) because their preparation requires the use of year-specific annual activity data such as the amount of fuel combusted or vehicle kilometers travelled in the inventory base year, and such activity data will only be available sometime after the end of the base year. Projected (or future-year) emissions inventories, on the other hand, are prospective (i.e., forecasts) because they are generated by “projecting” the emissions from a retrospective base inventory forward in time to a future year based on a set of assumptions about expected future changes to the economy, population and housing, on-road and off-road vehicle fleets, etc. along with information about any emissions changes that are expected to result from the phased implementation of existing emissions control legislation. When annual emissions levels are changing rapidly, the use of a retrospective inventory may not represent current emission levels as well as a projected inventory despite the greater uncertainties associated with a projected inventory due to the additional assumptions involved in its generation (e.g., Pan et al., 2014; Moran et al., 2020a). This was the reasoning used to choose each of the three projected future-year emission inventories listed above (in addition, importantly, to their availability): the RAQDPS023 became operational in 2021 at a time when the most recent available retrospective national emissions inventories were for base years 2015 or 2016 whereas overall North American emissions had been trending downwards for more than two decades (e.g., ECCC, 2022b; https://gispub.epa.gov/air/trendsreport/2021, last access: 18 March 2026).

After special emissions-processing-ready versions of the three projected annual national anthropogenic emissions inventories had been generated by ECCC for Canada and acquired for the U.S. and Mexico (U.S. EPA, 2019), they underwent extensive processing to prepare the SET4.0.0 emissions input files. As described in Sect. S2 the first step was to make some adjustments to the Canadian and U.S. inventories to account for additional and missing information. The inventories were then read by the SMOKE (Sparse Matrix Operation Kernel Estimation) emissions processing system (see https://www.cmascenter.org/smoke/, last access: 18 March 2026) to prepare the model-ready SET4.0.0 input emissions files. More details about this demanding and time-consuming task can be found in Sect. S2.1.

Table 9 compares the annual totals of model-ready anthropogenic emissions of the seven CAC species by country. These values were obtained from the full set of SET4.0.0 emissions files. It is clear from Table 9 that total U.S. emissions are larger than either Canadian or Mexican total emissions, but the relative size of the national emissions for the three countries depends on the pollutant species. For example, U.S. SO₂ emissions are twice as large as Canadian SO₂ emissions and U.S. NO_x emissions are four times larger than Canadian NO_x emissions, but U.S. emissions for the other five CAC species are all larger than the corresponding Canadian emissions by a factor of six or more. Interestingly, NO_x emissions from ocean shipping are comparable to land-based Canadian and Mexican NO_x emissions. Note too that this summary is based on SMOKE outputs and hence does not agree exactly with inventory totals (cf. Table 1 of the companion paper by Moran et al., 2026). One source of these differences is that the model domain does not contain all inventory emissions. This is especially true for Mexico, for which only 340 of Mexico's 2457 counties are completely or partially included in the model domain, but also applies to the U.S., for which 105 counties, including the state of Hawaii, and overseas territories such as Puerto Rico, the U.S. Virgin Islands, Guam, and American Samoa, are located outside of the model domain although their emissions are included in the U.S. NEI. Another source of these differences is the VOC definition: VOC in Table 9 is the sum of the emissions of the 11 model reactive VOC species (i.e., post-SMOKE) whereas VOC in Table 1 of Moran et al. (2026) corresponds to the bulk VOC emissions reported in the inventories, including unreactive species (i.e., pre-SMOKE). Third, the PM emission values in Table 9 are much lower than in the inventory totals in the other table because they have been reduced by scaling with the transportable-fraction (TF) factor to account for SGS removal processes after emission (Pace, 2005).

Table 9 Summary of annual, model-ready Canadian, U.S., and Mexican anthropogenic emissions and ocean shipping emissions (tonnes) for nominal 2021/22 period based on the SET4 emissions input files. Note that (i) inventory emissions from jurisdictions located outside of the model domain (e.g., southern Mexico, Hawaii, Puerto Rico) are not included in these totals, (ii) VOC emissions are calculated as the sum of the emissions of 11 model VOC species (excludes EOTH), (iii) fugitive PM_2.5 and PM₁₀ emissions have been reduced by transportable-fraction scaling but not by meteorological modulation, and (iv) in-flight aircraft, biogenic, biomass burning, lightning, volcanic, sea-salt, pollen and other biologicals, and aeolian dust emissions are not included.

Figure 3 shows the spatial distribution across the RAQDPS023 domain of mean July SET4.0.0 anthropogenic emissions of four chemical species: NO, VOC, NH₃, and PM_2.5. These spatial distributions are quite different due to variations in the relative contributions of different emissions source types for each species. For example, the contribution of ocean-going vessels to NO and PM_2.5 emissions is obvious in Figs. 3a and d. NO emissions are also associated with urban areas, highways, power plants, and agriculture. The spatial distribution of NH₃ emissions, on the other hand, is dominated by the agricultural sector, including fertilizer application for field crops and animal husbandry, which includes farm-level emissions but also concentrated animal feeding operations located in rural areas.

Figure 3 Spatial distribution of SET4.0.0 July anthropogenic emission of four chemical species: (a) NO (tonnes NO per grid cell), (b) VOC (tonnes VOC per grid cell), (c) NH₃ (tonnes NH₃ per grid cell), and (d) PM_2.5 (tonnes PM_2.5/grid cell). Note that “VOC” here is total model VOC, that is, the sum of the 11 emitted model VOC species (see Table 3), and PM_2.5 emissions are the sum of emissions of the six emitted PM_2.5 chemical components (SU, NI, AM, EC, POM, and CM).

3.11.2 Biomass burning emissions

Wildfires and other biomass burning (BB) sources such as prescribed burns can emit large quantities of atmospheric pollutants, including PM_2.5, NO_x, CO, VOC, and NH₃, and sometimes occur over many days in the case of long-lived fires (e.g., Akagi et al., 2011; Urbanski, 2014). In addition to their significant local impacts, smoke from wildfires can also be transported hundreds or even thousands of kilometres downwind (e.g., Jaffe et al., 2004; Pavlovic et al., 2016; Teakles et al., 2017). About 39 % of Canada and 34 % of the U.S. are covered by forest (https://data.worldbank.org/indicator/, last access: 18 March 2026), so wildfires started by lightning strikes or human activity are common in both countries (e.g., https://ciffc.ca/publications/canada-reports, last access: 18 March 2026; https://www.ncei.noaa.gov/access/monitoring/monthly-report/fire/, last access: 18 March 2026). Wildfires occur sporadically, however, and their initiation and evolution are difficult to predict. Moreover, wildfire emissions depend on many factors, including fire location, availability and condition of fuel, and past and current weather. Thus, in order to quantify and incorporate wildfire emissions into a CTM, wildfires occurring even in remote locations must be detected in near real time and their future hourly emissions predicted by the modelling system.

In 2013 ECCC, in collaboration with the Canadian Forest Service (CFS) and with contributions from the U.S. Forest Service (USFS), began to run an experimental version of the RAQDPS named FireWork, which was identical to the operational RAQDPS except for the addition of hourly BB emissions to the usual hourly anthropogenic emissions and the inline biogenic and sea-salt emissions (described in the next section). The BB emissions were calculated by combining NRT satellite-derived “hotspot” locations with estimates of fuel moisture, total fuel consumption, and fraction of smoldering vs. flaming combustion at these locations. This information was provided by the CFS Canadian Wildland Fire Information System (CWFIS), with hourly BB emission estimates produced using the Fire Emission Production Simulator (FEPS), a component of the USFS's BlueSky Modeling Framework (Lee et al., 2002; Anderson et al., 2004; Larkin et al., 2009; Pavlovic et al., 2016; Munoz-Alpizar et al., 2017). The CWFIS tracks current fire danger conditions and fire activity and characteristics across Canada and the U.S. (https://cwfis.cfs.nrcan.gc.ca/, last access: 18 March 2026), while FEPS manages data concerning fuel type, pollutant emission factors, and heat release characteristics of individual BB events (https://research.fs.usda.gov/pnw/products/dataandtools/fuel-and-fire-tools-fft, last access: 18 March 2026). Note that we refer here to BB emissions rather than wildfire emissions because of the possibility that large grass fires and prescribed burns may also be detected by satellite.

FireWork was then introduced operationally by ECCC in 2016 as a nearly identical copy of the RAQDPS, that is, identical except for the inclusion of NRT satellite-derived BB emissions (Munoz-Alpizar et al., 2017). Running these two very similar forecast systems allowed wildfire plumes to be identified by their PM_2.5 concentration (referred to as “fire-PM_2.5”), which were calculated by “brute force” as the difference between the FireWork and RAQDPS PM_2.5 fields. The hourly evolution of these smoke plumes could then be predicted and displayed (Pavlovic et al., 2016; https://weather.gc.ca/firework/index_e.html, last access: 18 March 2026). FireWork then underwent a number of upgrades after 2016 (see Table A1), including the implementation of the Canadian Forest Fire Emissions Prediction System (CFFEPS). CFFEPS incorporates the code for the CWFIS while introducing a more flexible and updateable alternative to the FEPS module along with a new fire plume-rise scheme (Chen et al., 2019a, b). Section S3 describes CFFEPS version 4.1 and the calculation of NRT BB emissions by the RAQDPS-FW023, the formal name of FireWork v023, which was implemented operationally at the same time as the RAQDPS023 (Chen and Menelaou, 2021).

3.11.3 Other natural emissions

Two other important sources of natural emissions are soil and vegetation (e.g., biogenic emissions) and the oceans (e.g., sea-salt emissions). Both of these sources of emissions were represented in the RAQDPS023 by emissions models.

Vegetation naturally releases numerous VOCs to the atmosphere as a byproduct of photosynthesis and other plant physiological processes (e.g., Fehsenfeld et al., 1992; Monson et al., 1995; Fuentes et al., 2000; Guenther et al., 2000). Many of these compounds, such as isoprene, monoterpenes, sesquiterpenes, and aldehydes and other oxygenated VOCs, are reactive and thus can influence atmospheric chemistry (e.g., Guenther et al., 2000; Duhl et al., 2008). Moreover, globally and in many regions, including eastern North America, these biogenic VOC emissions are larger in magnitude in total than anthropogenic VOC emissions (Roselle et al., 1991; Guenther et al., 1995, 2000; Pierce et al., 1998). Emissions of NO from microbial processes in certain soil types are a second type of biogenic emission (Williams et al., 1992a, b; Yienger and Levy, 1995; Stohl et al., 1996; Guenther et al., 2000). Soil NO emissions tend to occur in rural regions, where other anthropogenic sources of NO are small, and especially in agricultural areas where nitrogen fertilizer is applied since soil nitrogen availability is a controlling factor (Sha et al., 2021).

To calculate biogenic emissions at each chemistry time step the RAQDPS023 employed inline code based on a modified version of the Biogenic Emission Inventory System (BEIS) version 3.09, which was developed by the U.S. EPA (Vukovich and Pierce, 2002), together with a modified version of the Biogenic Emissions Landuse Database version 3 (BELD3), a North American vegetation database (Pierce et al., 2000; Zhang et al., 2018), and meteorological fields predicted by GEM. BEIS v3.09 calculates the emissions of four chemical species or groups of species – isoprene; monoterpenes; “other VOCs” (OVOC); and soil NO – using formulas from Geron et al. (1994). To save computer time, gridded summer and winter standard emissions fluxes (unit: g s⁻¹ grid cell⁻¹) were pre-calculated for standard mid-day meteorological conditions (ambient temperature of 30 °C, photosynthetic active radiation (PAR) photon flux density of 1000 µmol m⁻² s⁻¹) by summing vegetation-species-specific emission fluxes weighted by the vegetation-species fraction for each model grid cell (Geron et al., 1994; Guenther et al., 1994; Pierce et al., 1998). Winter emission fluxes for deciduous species were assumed to be zero. These standard biogenic emission flux fields were then treated like other time-independent geophysical fields and were input at the start of a RAQDPS023 simulation. A season-specific standard emission flux for each grid cell was then calculated during a simulation as the weighted average of the summer and winter standard emission fluxes using one of the following pairs of weights for the five phenological seasons of Zhang et al. (2002a): $\mathrm{1.0}/\mathrm{0.0}$, $\mathrm{0.66}/\mathrm{0.34}$, $\mathrm{0.25}/\mathrm{0.75}$, $\mathrm{0.0}/\mathrm{1.0}$, and $\mathrm{0.5}/\mathrm{0.5}$ for midsummer, autumn before harvest, autumn after harvest, winter, and spring seasons, respectively (cf. Sect. 3.9). The phenological season for each grid cell was determined based on the month and latitude (Brook et al., 1999a). GEM-predicted air temperature and solar angle fields were then used at each chemistry time step to calculate a light correction factor for isoprene and temperature correction factors for all three VOC species and soil NO. These correction factors accounted for the difference between the actual meteorological conditions and the standard meteorological conditions and were used to scale the season-specific standard biogenic emission fluxes. Lastly, emissions of the first two VOC species and soil NO were mapped directly to three ADOM-2 model species (ISOP, ALKE, and NO) while emissions of OVOC were split equally between ALKA and ALKE.

Note that three major updates were made to the standard BEIS v3.09 system described by Vukovich and Pierce (2002). First, the v3.09 biogenic emission factors have been replaced by newer BEIS v3.13 emission factors (Schwede et al., 2005). Second, standard emission fluxes for monoterpenes and OVOC were reduced for all pine species by factors of roughly five and two, and for all spruce species by factors of roughly five and three, as suggested by Stroud et al. (2010). And third, the original BELD3 vegetation database was updated over Canada to replace the three generic forest types (coniferous, deciduous, mixed) that had been specified for Canada with 30 BELD3 tree species based on the 2001 Canadian Forest Inventory (Beaudoin et al., 2014). One remaining limitation, however, that was addressed in later versions of BEIS was that soil NO emissions did not depend on soil moisture or precipitation (cf. Yienger and Levy, 1995).

The generation of sea spray aerosol over the world's oceans is an important PM source (e.g., De Leeuw et al., 2011). These sea-salt particles have a major influence on marine cloud properties and the radiation balance over the oceans, as well as coastal areas. A number of semi-empirical schemes are available for estimating sea-salt emission fluxes, but their estimates of these emissions at the global scale can vary by nearly two orders of magnitude (Textor et al., 2006; De Leeuw et al., 2011). One reason for this large uncertainty is the need for commonality in terms of assumed RH, reporting height, and particle size range.

Two schemes for sea-salt emission fluxes from the open ocean are available in MACH: the Smith et al. (1993) scheme, which is applicable to particles with wet diameters (RH = 80 %) in the 2 to 40 µm range, and the Gong-Monahan scheme (Gong, 2003), which is applicable to particles with wet diameters (RH = 80 %) in the 0.14 to 40 µm range. A third scheme for estimating sea-salt emission fluxes in the near-shore surf zone based on a single power law is also available in the model (De Leeuw et al., 2000, Eq. 4). This simple scheme is applicable to particles with wet diameters (RH = 80 %) in the 1.6 to 20 µm range. Note that sea-salt particles are assumed to equilibrate instantaneously with the ambient water vapour in all three schemes (Fitzgerald et al., 1998). As well, all three schemes are formulated in terms of 10 m wind speed and describe the dependence of particle number flux at 10 m height on wet particle size (RH = 80 %), which must then be used to calculate integrated dry particle mass flux estimates across each model size bin. Note that De Leeuw et al. (2011) suggest that sea-salt particle diameters at RH = 80 % are roughly double their diameters at RH = 0 %.

The RAQDPS023 used the latter two schemes, where the Gong-Monahan scheme was applied for any grid cell that includes at least some ocean while the surf zone scheme was used when the fraction of a grid cell covered by ocean was less than unity. For both schemes the resulting sea-salt emission flux estimates were scaled for each grid cell by the fraction of ocean coverage. Lastly, to calculate sea-salt emission fluxes for the two RAQDPS023 size bins, each bin was subdivided into 35 sub-bins in order to integrate the particle-size-dependent emission flux functions by particle radius.

Note, however, that some other sources of natural emissions, including lightning NO emissions, volcanic SO₂ emissions, pollen and other biological PM emissions, oceanic gas-phase emissions, and aeolian dust emissions, were not considered by the RAQDPS023. While all of these emissions affect atmospheric chemistry, their impacts on near-surface gas-phase and PM_2.5 concentrations over Canada are assumed to be small, especially when compared to biogenic, biomass burning, and sea-salt emissions, the three types of natural emissions that are considered. Lightning NO emissions occur intermittently above the Earth's surface. While they have little influence on surface concentrations locally, they can influence NO${}_{\mathrm{3}}^{-}$ wet deposition (e.g., Appel et al., 2011; Zhang et al., 2018b). Volcanic SO₂ emissions can affect atmospheric concentrations and deposition, but active volcanoes are not a significant source in populated areas of Canada or the U.S. although this is not true in Mexico (e.g., Fioletov et al., 2016). Pollen and other biological PM emissions are an important source of atmospheric PM and aeroallergens, but they are present mainly as particles larger than 2.5 µm in diameter, including pollen grains, which are typically larger than 20 µm, and fungal spores, which are typically found in the 1–10 µm range (e.g., Efstathiou et al., 2011; Sierra-Heredia et al., 2018; Subba et al., 2021). Oceanic emissions of gas-phase species such as DMS and isoprene contribute to background atmospheric chemistry, although DMS emissions are small at higher latitudes and marine isoprene emissions are small compared to terrestrial emissions (e.g., Bates et al., 1992; Khan et al., 2025). However, as discussed in Sect. 6.2 emissions of halogen species such as iodine can affect surface ozone levels. Aeolian dust emissions are discussed further in Sect. 3.12, but such emissions are very sporadic and are much less important in Canada than the U.S.

3.12 Meteorological modulation of fugitive PM emissions

Anthropogenic fugitive PM emissions are incidental soil dust emissions that originate from surface sources such as vehicles travelling on paved or unpaved roads and a number of industrial, construction, and agricultural activities like land-clearing, excavating, grading, demolition, tilling, and harvesting. The dominant chemical components of fugitive PM emissions are CM and POM (e.g., Boutzis et al., 2020). These emissions make up a major fraction of the anthropogenic primary PM emissions considered by regional AQ models. For example, in the 2020 Canadian APEI, fugitive dust emissions contributed 88 % of total primary PM_2.5 emissions and 94 % of total primary PM₁₀ emissions (ECCC, 2022a).

Unlike primary PM emissions from smokestacks and tailpipes (i.e., combustion processes), fugitive PM emissions can be directly affected by meteorological conditions since soil wetness or snow cover may reduce fugitive PM emissions. To account for such meteorological “modulation”, hourly inventory-based fugitive PM emissions were reduced in the RAQDPS023 at each time step during a model simulation according to predicted meteorological conditions in a surface grid cell if one of the following two meteorological thresholds was exceeded:

• i.

  grid-scale soil moisture (soil volumetric water content) greater than 10 %; or

• ii.

  grid-scale snow depth greater than 5 cm.

If either threshold was exceeded, then fugitive PM emissions from that grid cell for that time step were set to zero for all PM size bins and chemical components. This scheme was thus most active during rainy seasons and in winter. Note that in order to apply this scheme, input PM surface emissions had to have been separated and saved as non-fugitive and fugitive components when they were generated (see Sect. S2).

The choice of 10 % for the grid-scale soil wetness threshold was made based on sensitivity tests for several trial values followed by comparison of predicted near-surface CM concentrations with CM values measured in one Canadian and two U.S. PM_2.5 speciation networks (NAPS, CSN, IMPROVE). The U.S. operational AQ forecast system, the National Air Quality Forecast Capability (NAQFC), also applies a similar reduction of fugitive dust emissions by scaling PM emission fluxes by snow-cover fraction (Lee et al., 2017).

Note that this meteorological modulation scheme, which accounts for day-specific reductions of anthropogenic fugitive dust emissions by precipitation and snow cover, was independent of the land-cover-dependent TF scaling that was applied to anthropogenic fugitive dust emissions at the end of emissions processing (see Sect. S2). It should also be noted that natural aeolian (i.e., wind-blown) fugitive dust emissions, which were not considered by the RAQDPS023, can be important in arid or semi-arid parts of North America such as the desert southwest of the U.S. Such emissions are very episodic, however, and tend to occur in the spring and summer seasons under dry and very windy conditions (e.g., Orgill and Sehmel, 1976; Park et al., 2007; Foroutan et al., 2017; Kim et al., 2021). Nevertheless, these natural emissions may contribute over 50 % of PM_2.5 mass in the spring season in the southwestern U.S. and 20 %–30 % of PM_2.5 mass in the summer in the central and southeastern U.S. (e.g., Orgill and Sehmel, 1976; Hand et al., 2017; Tong et al., 2023), and the operational U.S. NAQFC does include an aeolian dust emissions scheme (Dong et al., 2016; Campbell et al., 2022).

3.13 Plume rise and vertical spread

Emissions from large smokestacks may disperse in the atmosphere quite differently from surface emissions. Such continuous pollutant plumes can rise well above the top of the emitting smokestack (plume rise) and then, depending on atmospheric vertical stability, display a number of types of dispersive behaviour, including coning, looping, fanning, lofting, and fumigating (e.g., Pasquill and Smith, 1983). Plume rise also exposes these pollutant plumes to different meteorological conditions at different elevations, including differing wind speeds and directions (referred to as vertical wind shear).

Pollutant point sources for which plume rise is modelled are often referred to as major point sources. In order to calculate plume rise, four major-point-source characteristics (known as stack parameters) are typically required: geometric stack height; geometric stack exit diameter; pollutant stack exit temperature; and pollutant stack exit flow speed (or stack volume flow rate). Major-point-source emissions were read by the RAQDPS023 from a separate file that also included stack-parameter values for each major point source (see Sect. S2). Information was also needed about atmospheric conditions at each stack location, including the vertical temperature profile, wind speed at the stack exit, PBL height, Monin-Obukhov length, and friction velocity. These values were all forecast by the GEM model (Table 2).

To calculate plume rise and vertical spread the RAQDPS023 employed a set of empirical formulations developed by Briggs (1984) that consider both buoyancy- and momentum-driven processes. Vertical spread represents the impact of SGS vertical diffusion associated with the mixing of the plume with its surrounding environment; it was expressed in the RAQDPS023 in terms of the vertical layers into which the emitted pollutant is injected. The implementation of the Briggs scheme in the RAQDPS023 has been described in detail by Akingunola et al. (2018). Note, though, that plume rise associated with wildfire plumes was treated differently (Sect. S3)

While the Briggs (1984) scheme provides estimates of mean injection height and plume vertical spread, one limitation of the treatment of major point source emissions by GEM-MACH (and other Eulerian CTMs) is that instantaneous horizontal and vertical mixing is assumed across the horizontal grid cell in which the point source is located and the vertical layers into which the pollutant is injected. This assumption usually has a minor impact on the far-field dispersion of the plume, but it does overestimate near-field plume diffusion and also changes plume chemistry by removing the distinction between SGS pollutant concentrations inside and outside of the plume. Note that a number of SGS plume-in-grid schemes have been proposed to address this limitation (e.g., Mathur et al., 1992; Kumar and Russell, 1996; Karamchandani et al., 2011), but the RAQDPS023 did not account for SGS plume chemistry.

4 Computational aspects

4.1 Time integration

For the time integration of a GDPS, RDPS, or RAQDPS simulation, time-stepping calculations start with the GEM dynamical core (Sect. S1.1). After the calculations for a dynamics time step have been completed, the GEM “physics” module (e.g., Mailhot et al., 1998; McTaggart-Cowan et al., 2019a) is then called to step through a set of parameterization schemes representing different atmospheric physical processes. Changes to GEM state variables that result from each physical process parameterization (i.e., tendencies) are treated by GEM as forcings that modify the dynamical equations using the sequential tendency-splitting technique for physics-dynamics coupling discussed by Gross et al. (2018). One exception, though, is the atmospheric radiative transfer scheme, which is computationally expensive and which, in both the RDPS 8.0.0 and RAQDPS023, was only called every fifth dynamic time step (i.e., every 1500 s).

In a GEM-MACH simulation the MACH module is called upon completion of the GEM physics module. In order for the RAQDPS023 to fit within its operational forecast “window” of 30 min wall-clock time, however, computational cost and execution time was reduced by using a chemistry integration time step of 900 s, three times longer than its dynamics integration time step of 300 s. This was necessary because even with this simplification the computational cost of MACH chemistry for one 900 s time step was approximately four times greater than the computational cost of GEM meteorology for three shorter 300 s time steps. When the start of an RAQDPS023 chemistry integration time step of 900 s coincided with the start of a dynamics integration time step, the MACH chemistry module was triggered after the GEM physics module had been executed. If the time steps did not align, though, the chemistry module was not called. This difference between dynamic and chemistry time steps was possible because chemistry processes were not permitted to impact dynamics/physics fields; that is, the RAQDPS023 was run with no feedback from chemistry fields to dynamics/physics fields, an approach referred to as one-way coupling (e.g., Grell et al., 2005). This is an acceptable approximation most of the time for tropospheric chemistry (though not for stratospheric chemistry). It should be noted, though, that the GEM-MACH code does support two-way coupling with feedbacks from chemistry fields to dynamics/physics fields as an option (e.g., Gong et al., 2015; Makar et al., 2015a, b).

As is common practice for CTMs, which must deal with a wide variation in time scales for different processes, process splitting is used in GEM-MACH to transform the coupled set of partial differential equations for tracer mass-conservation into a sequence of “smaller” problems where each process is described by a numerical operator (e.g., McRae et al., 1982; Marchuk, 1990; Pudykiewicz et al., 1997; Oran and Boris, 2000; Dimov et al., 2008; Wan et al., 2013). This operator-splitting approach simplifies the numerical solution of this system of equations and allows different specialized numerical solvers to be used for the chemistry time-step integration of different process operators. For example, as noted in Sect. 3.4 and 3.5, in order to integrate the gas-phase and aqueous-phase chemistry operators over a chemistry integration time step, each operator is treated using a specialized solver designed to handle a `stiff' set of ordinary differential equations, where much finer, variable “process-internal” time steps are employed. The solvers for each process are each applied once sequentially in a non-centered or non-symmetric manner, again to reduce computational cost (Dimov et al., 2008).

The operator sequence applied for the RAQDPS023 for each chemistry time step consisted of three groups of process operators: (i) emissions, vertical diffusion, and dry deposition; (ii) gas-phase chemistry; and (iii) aerosol processes. In the first process group, hourly anthropogenic and BB emission fluxes were input hourly and allocated to each chemistry time step while biogenic and sea-salt emission fluxes were calculated for each chemistry time step. Also in the first group, for each chemistry time step meteorological modulation was applied to fugitive PM emissions, plume rise was calculated for each point source and those elevated emissions were injected into the vertical column, dry deposition velocities were calculated, and then the GEM vertical diffusion scheme was applied with dry deposition fluxes used as bottom boundary conditions. Note that surface emission fluxes would normally be included in the specification of the bottom boundary fluxes, but as described at the end of Sect. 3.1, in the RAQDPS023 these emissions were instead injected into the two lowest model layers. The second operator group consisted of the gas-chemistry chemistry operator alone, which was applied for each chemistry time step but integration of this operator is performed using many smaller time steps in order to address numerical stiffness (Sect. 3.4). The third operator group consisted of many PM-related operators: one sub-group consisted of the operators for particle nucleation, condensation, and intersectional mass transfer of sulfuric acid and SOA precursor condensate; a second (cloud processing) sub-group included particle activation, gas-phase-aqueous-phase mass transfer and aqueous-phase chemistry, in-cloud scavenging and rainout, inorganic heterogeneous chemistry, and intersectional mass transfer; and the third sub-group consisted of particle coagulation followed by below-cloud precipitation scavenging of gases and particles by falling hydrometeors and hydrometeor evaporation.

In terms of computational cost, based on a number of timing tests for eight RAQDPS024 January and July hindcast cases (Moran et al., 2026) it was found that with chemistry turned off, GEM dynamics and physics each took about 45 % of run time, with model initialization accounting for most of the rest. Turning chemistry on (i.e., activating MACH), however, increased total run time by a factor of 4.4 on average. Interestingly, the computational burden of the set of chemistry operators, which were only called every third meteorological time step, increased total run time relative to GEM alone only by 107 %, that is, roughly a doubling in run time. However, with chemistry activated the cost of the dynamics calculations increased by a factor of 6.0 and the cost of the physics calculations, which include chemistry, increased by a factor of 3.4. The large increase in the cost of dynamics was due to the large increase in advection calculations due to the 58 extra chemical tracers that had to be advected every meteorological time step (Sect. 3.2, 3.4). Lastly, looking at the relative cost of the different chemistry operators we found that gas-phase chemistry operator took 23 % of the time required for chemistry and PM-related operators used 73 % of the time, with the calculation of plume rise, biogenic emissions, gas-phase dry deposition, and vertical diffusion accounting for the rest.

4.2 Initial conditions and boundary conditions

The RAQDPS023 followed the same initialization procedure used to provide meteorological initial conditions (ICs) for the RDPS 8.0.0, including the use of the same meteorological analysis files. Both systems employed a continuous data assimilation cycle and a four-dimensional incremental analysis update (4D-IAU) approach (Bloom et al., 1996; Buehner et al., 2015; Lei and Whitaker, 2016; CMC-RDPS-8.0.0, 2021a; Lu and Wang, 2021; Ritchie et al., 2022). In this approach, hourly analysis increments were applied from T−3 to T+3 h, where T=0 was the run start time, in order to maximize the dynamical balance between the mass and momentum fields through the application of relatively small but continuous forcing. The use of this scheme also allowed “recycling” of a number of unobserved model physics fields at T−3 h that were obtained from the previous operational RDPS 8.0.0 run, namely turbulence kinetic energy (TKE), turbulence regime, mixing length, friction velocity, PBL height, total condensate (combined liquid and ice) and cloud fraction from the Sundqvist scheme (Sect. S1.6), PBL cloud water and cloud fraction, and the flux enhancement factor from the MoisTKE scheme (Buehner et al., 2015; CMC-RDPS-6.0.0, 2018; CMC-RDPS-8.0.0, 2021b). Recycling these selected forecast physics fields reduced the spin-up time required for these fields to reach equilibrium in the new simulation. For example, cloud-related fields normally require a number of hours to reach equilibrium if initialized to zero, and the PBL-parameterization-related fields also contain information about PBL history that cannot be reconstructed diagnostically (Buehner et al., 2015).

Dynamical downscaling with one-way nesting was used to specify the meteorological LBCs (e.g., Fillion et al., 2010; Ritchie et al., 2022). Piloting files of hourly meteorological LBCs were provided from a RDPS 8.0.0 forecast to the RAQDPS023 forecast for the same start time, where the RAQDPS023 forecast followed the RDPS 8.0.0 forecast in wall-clock time, and the RAQDPS023 grid was a co-located subset of the RDPS 8.0.0 grid (Sect. 2.1). Hourly RDPS 8.0.0 forecast fields from the piloting files were then linearly interpolated in time to each 300 s GEM-MACH dynamics time step. The RAQDPS023 horizontal grid consisted of a large, interior “free” zone and a narrow (eight-grid-cell-wide) outer “blending” zone. A 12-grid-cell piloting zone also surrounded the RAQDPS023 free zone, so that it covered the blending zone but also extended another four grid cells outside of the RAQDPS023 grid domain. Following the approach of Thomas et al. (1998), open or inflow/outflow meteorological LBCs were applied at the RAQDPS023 horizontal grid boundary while a Davies-type boundary relaxation scheme (Davies, 1976) was applied across the blending zone after each dynamics time step to relax RAQDPS023 near-boundary values towards RDPS 8.0.0 piloting values using a cosine² cross-zone weighting function (see also Benoit et al., 1997). By smoothly blending RDPS 8.0.0 piloting values with interior values along the edges of the RAQDPS023 grid, the blending zone transmitted information about upstream atmospheric conditions external to the RAQDPS023 domain while also acting as a horizontal “sponge” zone to dampen spurious wave reflections from the lateral boundaries.

The meteorological upper boundary conditions, on the other hand, assumed a material surface with zero cross-boundary flow specified as a “truncated closed-top” boundary condition (Girard et al., 2014). Enhanced horizontal diffusion was also imposed at the top six model levels to reduce spurious wave reflections from the upper boundary by acting as a damping vertical “sponge layer” (CMC-RDPS-8.0.0, 2021b).

Chemical tracer ICs and LBCs are also required. The RAQDPS023 did not employ chemical data assimilation (CDA). It was instead run in a “perpetual forecast” mode, in which chemical IC fields for the start time of interest are copied from forecast chemical fields corresponding to that same time from the previous RAQDPS023 forecast launched 12 h earlier. Chemical LBCs (CLBCs), on the other hand, were specified using “climatological” seasonal vertical cross-sections of model species MMRs for each lateral boundary. These predetermined seasonal CLBCs were derived from a 2009 annual simulation of MOZART-4 (Model for OZone and Related chemical Tracers version 4), a global CTM that contains a detailed treatment of tropospheric inorganic chemistry and of some organic species (Emmons et al., 2010). The 2009 simulation was run on a 128 × 64 horizontal grid (approximately 2.8° × 2.8°) with 28 vertical levels reaching to 2 hPa. Global anthropogenic emissions of 20 model species were obtained from the ARCTAS inventory (Streets et al., 2006; D'Allura et al., 2011; http://bio.cgrer.uiowa.edu/arctas/emission.html, last access: 18 March 2026; https://www.acom.ucar.edu/gctm/mozart/subset, last access: 18 March 2026) and biogenic, biomass burning, marine, lightning, and volcanic emissions were also considered (Emmons et al., 2010). MOZART-4 does not have a complete stratospheric chemistry mechanism, so O₃, NO_x, and CO mixing ratios were constrained to climatological values from 50 to 2 hPa, and by relaxing to the climatology with a time scale of 10 d from 50 hPa to the tropopause. An annual set of six-hourly chemistry outputs from MOZART-4 was mapped to GEM-MACH model species, divided and averaged over each season, and horizontally and vertically interpolated to the RAQDPS023 grid to obtain four sets of seasonal CLBCs (Pendlebury et al., 2018). While this approach for providing CLBCs was an improvement over schemes used in earlier versions of the RAQDPS, where constant values or static vertical profiles were used at the lateral boundaries, it was still computationally simple and efficient.

4.3 Tracer advection, shape preservation, and mass conservation

RAQDPS023 prognostic chemical tracers such as O₃, NH₃, PM_2.5-EC, and PM_cf-POM were treated by the GEM dynamical core in the same way as the two RDPS 8.0.0 and RAQDPS023 prognostic meteorological tracers (water vapour, cloud condensate) described in Sect. S1.1. The GEM semi-Lagrangian advection scheme was thus also used to advect 58 GEM-MACH chemical tracers: 42 gas-phase tracers (Table 3) and 16 (= 2 × 8) PM size bin-chemical component tracers (Sect. 3.2), where the same two-step ILMC monotonicity correction used by GEM for meteorological tracers to impose shape preservation and avoid the creation of artificial local minima or maxima by the advection scheme (cf. Sect. S1.1) was also used for the chemical tracers. In addition, the Bermejo and Conde (2002) global mass fixer with LAM boundary-flux estimation (Diamantakis and Flemming, 2014; Aranami et al., 2015) that was applied to the meteorological tracers in the RDPS 8.0.0 and RAQDPS023 to impose mass conservation (Sect. S1.1) was also applied in the RAQDPS023 to the chemical tracers (de Grandpré et al., 2016).

One concern with global mass adjustments made to tracer species fields to conserve mass is that they have the potential to introduce artificial numerical diffusion by being imposed at locations far from locations where mass-conservation violations occurred. Diamantakis and Flemming (2014) found, however, that global mass tracer conservation schemes could improve forecast accuracy, and de Grandpré et al. (2016) recommended the use of the ILMC scheme for shape preservation in part due to its good mass-conservation properties, which then result in smaller mass adjustments for mass conservation. Note that one additional complication with the advection of a suite of chemical tracers is that they are often chemically coupled, so that the imposition of shape preservation and mass conservation on individual chemical species as was done in the RAQDPS023 may result in the generation of chemical imbalances (e.g., Lauritzen and Thuburn, 2012; Diamantakis, 2013), which then require equilibration by the gas-phase chemistry operator (Sect. 3.4).

4.4 Computer hardware and code parallelization

When the RAQDPS023 was implemented operationally in late 2021, it was run on the ECCC high-performance computing (HPC) system then available. The heart of this system consisted of two Cray XC50 massively parallel supercomputers, each with 1330 shared-memory compute nodes and 53 200 compute cores, where each compute node had 192 GB of random-access memory. Having two identical computers available for operational use provides redundancy and allows hardware maintenance, software upgrades, and various system diagnostics and tests to be performed on one computer while operational forecasts continue to be run on the other computer. As well the two supercomputers were located in separate halls to reduce physical risk. In June 2022 the two Cray supercomputers were replaced by two new IBM ThinkSystem SD650v2 direct-water-cooling supercomputers. Each of these new machines has 1496 shared-memory compute nodes and 119 680 compute cores, where each compute node has 512 GB of random-access memory. The RAQDPS023 was migrated without any algorithmic upgrades to this new HPC system and was re-named RAQDPS024 in June 2022 (see Table A1).

The GEM and GEM-MACH Fortran2003 codes support both the Message Passing Interface (MPI) and Open MultiProcessing (OpenMP) parallelization paradigms. This allows the RAQDPS023 to employ hybrid parallelization, that is, distributed-memory and shared-memory parallel programming techniques, respectively, to take advantage of the large number of available compute nodes and cores. Domain decomposition has been used to implement MPI parallelization, where the RAQDPS023 horizontal grid was divided into 30 subgrids or “tiles” in the x-direction (longitude) and 36 tiles in the y-direction (latitude) for a total of 1080 tiles. Integration was then performed simultaneously on these tiles on 1080 compute nodes, where for each tile, calculations were carried out sequentially across each tile or subgrid from south to north on a set of XZ vertical “slices”. Two OpenMP threads were used to speed up the integration.

4.5 Forecast run strategy

The RAQDPS023 was run two times per day for 72 h simulations starting (nominally) at 00:00 and 12:00 UTC. These runs, which had to follow the 00:00 and 12:00 UTC RDPS 8.0.0 (and GDPS 8.0.0) runs that were needed to provide piloting fields (Fig. 1), started at 03:00 and 15:00 UTC wall-clock time. Each RAQDPS023 run, including pre- and post-processing tasks, took approximately 90 min of wall-clock time, including approximately 30 min for the 75 h GEM-MACH simulation; the extra three hours included the $T=-\mathrm{3}$ h start for the 4D-IAU initialization procedure (Sect. 4.2). Forecast results were available for the 00:00 and 12:00 UTC runs at approximately 04:30 and 16:30 UTC, respectively, or close to local midnight and local noon in both Atlantic Standard Time and Eastern Standard Time. The independent runs of the RAQDPS-FW023 launched and finished at similar times, with the one extra complication that RDPS 8.0.0 72 h forecasts were also needed by CFFEPS v4.1 (which must be run before the RAQDPS-FW023) for estimating BB emissions for the next 72 h (see Sect. S3). Note that the start times and run times for the RAQDPS024 and RAQDPS-FW024 systems were very similar to those for the RAQDPS023 and RAQDPS-FW023.

5 Outputs, post-processing, downstream systems, and performance evaluation

5.1 Outputs

Routine RAQDPS023 (and RAQDPS-FW023) hourly outputs included a large number of surface chemical fields, among them 41 gas-phase species MMR fields (first 41 species in Table 3), eight PM_2.5 and eight PM_cf chemical component MMR fields, PM_2.5 and PM₁₀ total mass concentration fields, 21 dry deposition flux fields (SO₂, H₂SO₄, NO, NO₂, O₃, HNO₃, HONO, PAN, RNO₃, NH₃, H₂O₂, HCHO, NO_y, SU, NI, AM, CM, EC, POM, SOM, SS, where SU=SU1+SU2, etc.), 12 wet deposition flux fields (HSO${}_{\mathrm{3}}^{-}$, SO${}_{\mathrm{4}}^{=}$, NO${}_{\mathrm{3}}^{-}$, EC, POM, SOM, CM, SS, NH${}_{\mathrm{4}}^{+}$, H₂O₂, ROOH, H₂O), and four biogenic emission flux fields (ISOP, ALKA, ALKE, NO). RAQDPS023 forecasts of NO₂, O₃, and PM_2.5 surface fields in particular were used as guidance by operational AQ forecasters at ECCC and elsewhere, were provided as inputs to several post-processing systems, and were made available to the public through ECCC's public weather website (see https://weather.gc.ca/firework/index_e.html, last access: 18 March 2026). It should be noted that the PM_2.5 and PM₁₀ total mass values were calculated as the sum of seven chemical components (SU, NI, AM, EC, POM, SOM, CM); they thus correspond to dry PM_2.5 and dry PM₁₀ total mass (i.e., without aerosol water) without a contribution from sea salt, thus emphasizing values over land and without a dependence on RH.

For public outreach on the weather website, ECCC also provides forecasts of the AQHI as a key product to communicate the total health risk of a mixture of air pollutants and also to disseminate warnings during forecasted periods of poor air quality (https://weather.gc.ca/mainmenu/airquality_menu_e.html, last access: 18 March 2026). The AQHI is a multi-pollutant, additive, no-threshold, health-based, hourly AQ index with a range (in whole numbers) from 0 to 10+. It was developed from daily time-series analysis of air pollutant concentrations and mortality data and is calculated as a weighted sum of NO₂, O₃, and PM_2.5 surface concentrations (Stieb et al., 2008; Trieu et al., 2020). Two interpretative health messages are provided with a 72 h AQHI forecast, one for the general population and one for at-risk groups, including children, the elderly, and individuals with heart or breathing problems. While O₃ and PM_2.5 are known to have direct health effects (e.g., Murray et al., 2020), NO₂ is included in the index largely as a proxy for other, non-measured pollutants with health impacts that are emitted by combustion processes (Stieb et al., 2008). It is worth noting that AQHI advisories and warnings in Canada in the warm season are increasingly being driven by elevated PM_2.5 levels associated with wildfire smoke forecasts.

5.2 Post-processing tasks and downstream systems

Figure 1 also shows the data flow between the RAQDPS023 and two operational downstream post-processing systems, the Updateable Model Output Statistics for Air Quality (UMOS-AQ) system and version 2.0.0 of the Regional Deterministic Air Quality Analysis (RDAQA) system. UMOS-AQ is a statistical AQ post-processing package for bias correction that can compensate for systematic model errors, account for unresolved SGS phenomena, and be updated frequently. An operational UMOS post-processing package for meteorology has been used by EC since 1995 to forecast meteorological predictands such as surface temperature and probability of precipitation (Wilson and Vallée, 2002, 2003). UMOS-AQ uses the UMOS framework to combine multiple sources of AQ-related information: AQ forecasts; current AQ measurements; meteorological forecasts; and physical variables (e.g., solar flux, sine of scaled Julian day). AQ forecasts made by the RAQDPS023 provided the operational UMOS-AQ package with point-specific hourly surface forecasts of NO₂, O₃, and PM_2.5 concentrations at Canadian AQ measurement station locations. The UMOS-AQ package itself consists of a large set of multi-variate linear regression (MLR) equations, where one MLR equation is generated per Canadian AQ station per pollutant per season per forecast hour per 00:00/12:00 UTC run (Antonopoulos et al., 2011; Moran et al., 2014). These equations are regenerated every week using the most recent four months of model AQ forecasts and station measurements, which permits ongoing adjustments to account for changes in model forecast skill due to many causes, such as forecast system upgrades, systematic seasonal variations, and unexpected emissions changes such as accompanied the COVID-19 pandemic in 2020 (e.g., Griffin et al., 2020; Mashayekhi et al., 2021). UMOS-AQ-derived values of NO₂, O₃, and PM_2.5 are used to calculate the location-specific AQHI forecast values that are disseminated to the Canadian public (see https://weather.gc.ca/mainmenu/airquality_menu_e.html, last access: 18 March 2026).

A second operational AQ post-processing system, the RDAQA system, provides objective analyses (OA) of hourly North American surface concentration fields for O₃, NO₂, NO, PM_2.5, PM₁₀, SO₂, and the AQHI (Robichaud and Ménard, 2014; Robichaud et al., 2016; Robichaud, 2017; Ménard, 2021). These surface OA fields are produced by combining recent hourly AQ surface measurements at station locations with hourly gridded model surface forecast fields in a kind of measurement-model fusion (e.g., Fu et al., 2022). AQ measurements are in general accurate and unbiased, but they have limited spatial coverage and may not always be spatially representative, especially in locations with complex topography or heterogeneous land cover. Model forecast fields, on the other hand, may have inaccuracies and biases, but they also provide complete spatial coverage and areal representativeness and should account for the meteorology, physics, and chemistry of air pollution. An objective analysis that is an optimal combination (known as optimal interpolation or OI) of both of these information sources will lead to a significant improvement of the coverage and accuracy of analyzed air pollution spatial patterns (Robichaud et al., 2016). Hourly analyses are prepared because many pollutants, including O₃ and NO₂, have marked diurnal variations and also respond to synoptic meteorological variations.

A new version of the RDAQA post-processing system, version 2.0.0, was introduced operationally at ECCC in December 2021 at the same time as the RAQDPS023 (Ménard, 2021). RDAQA 2.0.0 employed a new formulation of the error covariances with a slightly modified analysis solver. The error covariances were derived from ensembles of RAQDPS023 model simulations for each hour of the day and captured non-homogeneous, non-isotropic error correlations related to terrain-dependent features such as mountains, valleys, water surfaces, and emissions (Ménard et al., 2020; Ménard, 2021). These simulation ensembles were “climatological” and were composed of about 60 RAQDPS023 forecast runs from a previous two-month development period. It was assumed that error correlations for these ensembles would be sufficient to capture stationary and terrain-dependent structures associated with terrain features and surface emissions sources. RAQDPS023 hourly gridded forecast fields were used by the RDAQA 2.0.0 as first-guess fields. The resulting linear combinations of hourly AQ measurements and RAQDPS023 forecasts minimized the analysis (or combination) error variance. Figure 4 shows an example of hourly output fields from the RDAQA 2.0.0 for PM_2.5 at 16:00 UTC as a four-panel analysis.

Figure 4 Four-panel presentation of RDAQA 2.0.0 results for PM_2.5 surface concentration (µg m⁻³) for 28 September 2023 at 16:00 UTC: (a) RAQDPS023 PM_2.5 forecast; (b) RDAQA 2.0.0 PM_2.5 objective analysis; (c) PM_2.5 analysis increment (correction to the model forecast); and (d) PM_2.5 surface observations.

The availability of hourly deposition flux forecast fields from the RAQDPS023 also allowed seasonal and annual dry and wet deposition fields to be accumulated on an ongoing basis. These fields could then be used to examine deposition trends and ecosystem impacts (e.g., Vet and Ro, 2008; Makar et al., 2018a, 2025). These gridded multi-month deposition fields can also be used by the ECCC ADAGIO (Atmospheric Deposition Analysis Generated by optimal Interpolation from Observations) tool, which is being developed to produce OAs of dry, wet, and total deposition of acidifying pollutants over North America (Robichaud et al., 2020) in parallel to some similar U.S. efforts (Schwede et al., 2019; Fu et al., 2022).

Wildfire smoke plume forecasts were also generated as a post-processing step. As described in Sect. 3.11.2 the RAQDPS023 served as the base platform for the ECCC RAQDPS-FW023 (“FireWork”) wildfire forecast system, which augmented the air pollutant emissions considered by the RAQDPS023 with NRT emissions from wildfires and other large BB events (Pavlovic et al., 2016; Munoz-Alpizar et al., 2017; Chen et al., 2019a; Chen and Menelaou, 2021). Both RAQDPS023 versions were run twice daily at 00:00 and 12:00 UTC for 72 h periods. Hourly forecasts of the transport and evolution of wildfire smoke plumes were then calculated as the difference between the RAQDPS-FW023 and RAQDPS023 PM_2.5 forecast fields. Four post-processing products for wildfire smoke plumes were prepared and posted to an ECCC public website (https://weather.gc.ca/firework/index_e.html, last access: 18 March 2026): (i) 72 h animations of hourly ground-level smoke concentration fields (µg m⁻³); (ii) 72 h animations of hourly, vertically-integrated total column smoke fields (kg m⁻²) (this product was discontinued on 11 June 2024); (iii) 24 h mean ground-level smoke concentration fields for Days 1, 2, and 3; and (iv) 24 h maximum ground-level smoke concentration fields for Days 1, 2, and 3. Figure 5 shows examples of two of these products. RAQDPS-FW023 forecasts were also made available on the WMO North American Regional Vegetation Fire and Smoke Pollution Warning and Advisory Centre (RVFSP-WAC) website, which is operated by ECCC (https://hpfx.collab.science.gc.ca/~svfs000/na-vfsp-was/public/dist/home/introduction, last access: 18 March 2026).

Figure 5 RAQDPS-FW023 12:00 UTC forecast on 28 June 2023 of contribution of smoke from wildfires to mean surface PM_2.5 concentrations (µg m⁻³) for hours (a) 00:00–24:00 and (b) 24:00–48:00.

5.3 Performance evaluation

The operational performance of the RAQDPS023 was routinely monitored and reported internally (e.g., Moran et al., 2021a), but only the small number of pollutants whose NRT measurements were routinely available from surface networks could be considered. These pollutants include the three AQHI component pollutants, NO₂, O₃, and PM_2.5 total mass, plus NO, SO₂, and PM₁₀ total mass. These same six pollutants were considered by the RDAQA 2.0.0 system described above. Note that performance evaluation statistics are usually based on point measurements, but analysis increments from the RDAQA 2.0.0 can be viewed as regional errors (e.g., Fig. 4). The companion paper by Moran et al. (2026) has examined RAQDPS023 performance for five years: the first year of RAQDPS023 operational forecasts and four years of retrospective annual simulations for the 2013–2016 period. For the 2021/22 forecasts of NO₂, O₃, and PM_2.5 total mass made by the RAQDPS023, the seasonal and annual evaluation scores for O₃ hourly forecasts were generally the highest, followed by NO₂ scores and then PM_2.5 scores. Another important finding was that predicted monthly mean PM_2.5 total mass was biased low in all months in 2021/22.

Other pollutants are also monitored, but these measurements are only available after a lag of months or even years to allow time for laboratory analysis and quality assurance/quality control (QA$/$QC) measures and hence can only be used for retrospective performance evaluations. The companion paper by Moran et al. (2026) also examined RAQDPS023 performance using this larger suite of species measurements, including other gas-phase species, PM_2.5 chemical components, and concentrations in precipitation, for the 2013–2016 retrospective simulations, which used year-specific anthropogenic emissions in place of the SET4.0.0 emissions set. Access to the much more comprehensive measurement data set provided additional insights into model performance. One insight was that one PM_2.5 chemical component (SU) was consistently underpredicted for the four years while two others (EC, SS) were consistently overpredicted. O₃ predictions in the spring were also biased low while isoprene predictions were biased high. These findings all suggest possible directions for model improvement. For more details see the companion paper by Moran et al. (2026).

ECCC has also been working with other agencies since 2017 to exchange, evaluate, and compare operational AQ forecasts for North America. ECCC has built an automated verification system to receive, ingest and compare these daily forecasts against surface measurements (Pavlovic et al., 2018; Moran et al., 2020b). Three regional AQ forecast models and three global AQ forecast models from five agencies are currently participating in this initiative. A performance evaluation report containing North American monthly evaluation statistics for these six models is produced every quarter by ECCC and is then made available through the WMO Global Air Quality Forecasting and Information System (GAFIS) initiative website (https://community.wmo.int/en/activity-areas/gaw/science-for-services/gafis, last access: 18 March 2026) on webpage https://hpfx.collab.science.gc.ca/~svfs000/na-aq-mm-fe/dist/ (last access: 18 March 2026). This effort is discussed in more detail in Moran et al. (2026).

6 Further development and future work

The previous sections have described the RAQDPS023 and RAQDPS-FW023 operational AQ forecast systems that became operational on 1 December 2021 (Table A1). Moran et al. (2026), the companion paper to this paper, documents both strengths and weaknesses of these two systems, including some model weaknesses that future system versions might address. During the preparation of these two papers, however, a major new version of the RAQDPS, RAQDPS025, which does at least partially address some of these weaknesses, including PM_2.5 total mass underpredictions, was implemented operationally on 11 June 2024 (Table A1). This section briefly summarizes the main updates and innovations introduced in the RAQDPS025 and then outlines some possible improvements for future system versions.

6.1 RAQDPS025 updates and innovations

The RAQPS025 system is documented in CMC-RAQPS025 (2024). As described below it includes two innovations related to meteorology, eight innovations related to chemistry, and one innovation related to input emissions. So-called “unit” tests were performed during the development of this new version to quantify the impact of each innovation, similar to the approach described by Foley et al. (2010). Five of these innovations were judged to have a “high” impact, five were judged to have a “medium” impact, and one was judged to have a “low” impact (CMC-RAQPS025, 2024).

The key acceptance criterion for the adoption of a new operational forecast system at ECCC is that it should improve, or at least not degrade, overall system performance regardless of whether the proposed system is believed to deliver improvements to science representations or numerical schemes or system architecture. An ECCC internal management committee is responsible for approving parallel and operational implementations of new forecast system versions. If preliminary approval is given by the committee based on evaluation results from several retrospective evaluation periods, the new forecast system is run side-by-side (i.e., in parallel) with the existing forecast system for a minimum of two months before results of evaluations for the parallel period are presented to the committee and final approval for implementation is either given or withheld. The RAQDPS025 was run in parallel with the RAQDPS024 from 13 February to 10 June 2024.

6.1.1 Updates to GEM and GEM-MACH computer codes

The RAQDPS023 was based on the GEM 5.1.0 and GEM-MACH v3.1.0.0 computer codes, and the RAQDPS024 was based on the GEM 5.1.2 and GEM-MACH v3.1.1.2 computer codes. The RAQDPS025 is based on the GEM 5.2.1 and GEM-MACH v3.2.1.1 computer codes (Table A1). Many of the following innovations were introduced with GEM-MACH v3.2.1.1.

6.1.2 Updates to meteorological piloting model

One ongoing challenge for operational forecasting at the Canadian Meteorological Centre (CMC) is providing the computer and human resources required (i) to maintain the large number of complex weather and environmental forecast systems that are being run operationally and (ii)  to implement a suite of system updates at regular intervals while addressing the complicated interdependencies of these systems. As noted in Sect. 2 the regional 10 km RAQDPS023 was dependent on both the global 15 km GDPS 8.0.0 and the regional 10 km RDPS 8.0.0. In June 2024 a new version of the GDPS, GDPS 9.0.0, with horizontal grid spacing reduced to 10 km was implemented operationally, and at the same time the RDPS 8.0.0 was retired as an independent forecast system without a replacement (CMC-RDPS-9.0.0, 2024). As a consequence it was necessary to replace the RDPS 8.0.0 as the meteorological piloting model for the RAQDPS025 with the GDPS 9.0.0, but there was no loss of horizontal resolution. In principle this change also has the added benefit of allowing access to the GDPS 9.0.0 ozone field, which is updated by satellite data assimilation.

6.1.3 Merger of RAQDPS and RAQDPS-FW forecast systems

In order to identify and track BB smoke plumes by means of the “brute-force” method (Sect. 3.11.2) the ECCC AQ forecast system architecture required the same modelling system to be run twice at each forecast time, once with BB emissions (RAQDPS-FW023) and once without BB emissions (RAQDPS023). However, other approaches to track BB smoke plumes with only a single run are possible. One dynamic approach is “tagging”, in which BB emissions are tracked separately during a model simulation as “tagged” species (e.g., Wagstrom et al., 2008; Wang et al., 2009; Samaali et al., 2011; Grewe, 2013). Another approach is to include BB emissions but identify the presence of BB smoke at the post-processing stage by using known chemical characteristics of BB smoke, including both speciated PM components and gas-phase species (e.g., Akagi et al., 2011; Urbanski, 2014; Andreae, 2019; Hayden et al., 2022), to identify the locations of ambient concentrations of multiple species predicted by the model that match those chemical characteristics of BB smoke. Either approach would allow the existing system architecture to be streamlined by removing the need for two simulations with two versions of the same forecast system. However, two issues that must be considered are the different “signatures” of BB emissions associated with different ecosystems and the formation of secondary air pollutants such as O₃ and SOA in BB smoke plumes in addition to the primary pollutants emitted by BB.

The second approach was adopted in the RAQDPS025 to track BB smoke plumes in a single RAQDPS025 simulation, where the ranges of two ratios of three primary PM_2.5 components (EC, POM, CM) were used together to identify BB-dominated smoke plumes (CMC-RAQPS025, 2024). The lower and upper thresholds of these two ratios that identify a BB-dominated PM_2.5 concentration were determined from an analysis of chemical speciation profiles of PM emissions associated with different types of wood combustion, including wildfires. This tracer-ratio approach does represent a change in strategy since the previous brute-force approach yielded PM_2.5 concentration values due only to the contribution of BB emissions (i.e., fire-PM_2.5) whereas the PM_2.5 concentration values produced by the new approach also include the contributions of emissions from other emission source types. This loss of source differentiation will be less of an issue, however, when BB emissions are the dominant PM_2.5 source. Note that in contrast to the RAQDPS023 and RAQDPS-FW023, the RAQDPS025 is referred to at ECCC as a unified AQ forecast system.

6.1.4 Updates to representation of vertical diffusion in the stable planetary boundary layer

The design philosophy for the development of the RAQDPS has been to employ GEM parameterizations without modification whenever possible. As noted in Sect. 3.1, however, the RAQDPS023 treated the vertical diffusion of chemical tracers in stable PBLs slightly differently than the vertical diffusion of meteorological tracers. GEM did impose two limiting values related to its turbulence parameterizations: a minimum Obukhov length (Sect. S1.3) and a maximum eddy dissipation length scale (Sect. S1.4). Many CTMs, on the other hand, also specify a minimum threshold value for the vertical eddy diffusivity of pollutant tracers (K_z) to avoid overpredictions of surface concentrations under stable conditions (e.g., McNider and Pour-Biazar, 2020), and a number of studies have examined CTM forecast sensitivity to the choice of this minimum value (e.g., Pleim, 2007; Jin et al., 2010; Makar et al., 2014; Li and Rappenglueck, 2018; Jia and Zhang, 2021; Jiang et al., 2024). The RAQDPS023 employed a minimum K_z value of 0.1 m² s⁻¹. Another CTM, the CMAQ (Community Multiscale Air Quality) model extends this approach by considering land-surface-type-dependent minimum vertical eddy diffusivity values, which are specified to be higher over urban areas (Li and Rappenglueck, 2018). This approach is appealing since urban areas have large surface roughness values and are typically also areas of high surface emissions in which surface concentrations are very sensitive to the exact balance between near-ground turbulent vertical mixing and surface emissions (Makar et al., 2014; Li and Rappenglueck, 2018). However, the specification of the allowable range of minimum eddy diffusivity values in the CMAQ scheme is still arbitrary.

Makar et al. (2021) tested a novel vertical diffusion parameterization in which the near-surface vertical diffusion of on-road vehicle emissions is modulated by vehicle-induced turbulence (VIT), which is calculated from traffic density and vehicle speeds. Like the CMAQ scheme this VIT scheme can increase vertical mixing in urban areas for stable conditions, but it also offers a physical basis for this increase and it will affect non-urban areas with high levels of traffic as well (e.g., rural highways). The VIT scheme has been implemented in the RAQDPS025. As might be expected it has a greater impact in the winter months when stable conditions are more common. A co-benefit of adding this new scheme was the elimination of the arbitrary lower limit on PBL height of 100 m and the imposed injection of area emissions into the two lowest model layers instead of the first layer alone that was used in the RAQDPS023 (Sect. 3.1).

6.1.5 Updates to chemical lateral boundary conditions

CLBCs are important for regional AQ models and may even be controlling for those species whose lifetimes are greater than a few days, which is often the time needed for an air parcel to traverse a regional model domain (e.g., Brost, 1988). Among these longer-lived species are (i) tropospheric O₃, which has a lifetime of a week to several months and a globally averaged tropospheric lifetime of 22–23 d (e.g., Archibald et al., 2020), (ii) PM_2.5, which has a tropospheric lifetime on the order of days to a week (e.g., Heintzenberg, 1989; Pandis et al., 1995), and (iii) CO, whose lifetime is 30–90 d (e.g., Seinfeld and Pandis, 2016). The western, northern, and eastern lateral boundaries of the RAQDPS023 continental domain are located in low-emissions regions over the Pacific Ocean, Arctic Ocean, and Atlantic Ocean, respectively, while a portion of the southern boundary is located over Mexico roughly 2500 km from the Canada-U.S. border (Fig. 2). This combination of low emissions and distance will help to reduce the contribution of CLBCs to elevated pollutant levels over Canada. On the other hand, O₃ and CO have high background levels, and those levels can vary significantly with time at the lateral boundaries. For O₃ those variations may be due to processes such as stratospheric O₃ intrusions and trans-Pacific transport of Siberian wildfire plumes (e.g., Heald et al., 2003; Flemming et al., 2017; Pendlebury et al., 2018). PM_2.5 levels can also vary with time over the oceans due to trans-oceanic transport of wildfire plumes from Siberia and desert dust from the Sahara and from northern China, Mongolia, and Central Asia (e.g., McKendry et al., 2008; Teakles et al., 2017; Tang et al., 2021).

As described in Sect. 4.2, the RAQDPS023 used “climatological” seasonal CLBCs based on results from a 2009 annual simulation with the older MOZART4 offline global CTM. For the RAQDPS025 new climatological monthly CLBC files were generated for the ADOM-2 gas-phase species by calculating monthly averages from outputs for the 2012-2019 period extracted from a 20-year global simulation made with version 6 of the Community Atmosphere Model with chemistry (CAM-chem) (Emmons et al., 2020; Tilmes et al., 2022; CMC-RAQPS025, 2024). Year-specific anthropogenic emissions for these annual simulations were obtained from CAMS-GLOB-ANT v5.1, a recent multi-year global emissions inventory (Soulie et al., 2023, 2024), and natural emissions from biogenic, biomass burning, marine, lightning, dust, and volcanic sources were also considered (Emmons et al., 2020).

While such climatological CLBCs do not address the episodic nature of the large-scale distributions of longer-lived pollutants such as O₃, CO, and PM_2.5, they can still be broadly representative of typical chemical inflows. The alternative is to supply time- and space-varying CLBCs based on forecasts from a compatible global chemical weather model. Such a system is currently under development at ECCC, but it is not yet ready for operational deployment.

6.1.6 Updates to representation of below-cloud scavenging of size-resolved PM

As described in Sect. 3.10 the treatment of below-cloud scavenging of PM in the RAQDPS023 was based on Slinn (1984). More recent studies by Wang et al. (2014a, b) have proposed an updated and unified parameterization of below-cloud wet removal of size-resolved PM by different precipitation types in all seasons. This newer parameterization was tested in a research version of GEM-MACH (Ghahreman et al., 2024) and was then implemented operationally in the RAQDPS025.

As well the RAQDPS023 used a fixed threshold temperature of 0 °C to define the transition between liquid cloud droplet scavenging and solid snow scavenging. In the RAQDPS025 the phase state for cloud droplets predicted by GEM cloud physics, which includes mixed-phase as well as liquid and solid hydrometeors, has been used instead (Ghahreman et al., 2024). This better aligns the treatment of below-cloud scavenging of pollutants in GEM-MACH with the treatment of cloud processes in GEM.

6.1.7 Updates to representation of gas-phase dry deposition

Three minor updates were implemented to the gas-phase dry deposition code in the RAQDPS025 (CMC-RAQPS025, 2024): (i) a correction was made to the value of the Schmidt number (changed from 0.84 to 0.62); (ii) the algorithm for calculating the ratio of the molecular diffusion coefficient for water vapour to that of a depositing pollutant species was updated; and (iii) a correction was made to the calculation of aerodynamic resistance R_a. In addition, the monthly climatology of satellite-derived, gridded LAI described by Zhang et al. (2021), which provides high-resolution (500 m), latitude-longitude-dependent descriptions of variations in this important vegetation phenological attribute, was adopted in the calculation of the surface resistance term R_c. Note that the use of monthly LAI fields reduces but does not remove the dependence on five phenological seasons that was described in Sect. 3.9.

6.1.8 Inclusion of sea salt mass in PM_2.5 total dry mass calculation

As noted in Sect. 5.1 the RAQDPS023 calculated PM_2.5 total mass as the sum of seven chemical components (SU, NI, AM, EC, POM, SOM, CM). This quantity corresponds to dry PM_2.5 total mass (i.e., without aerosol water) without a contribution from sea salt. The SS contribution had been excluded in earlier versions of the RAQDPS due to large overpredictions that had been observed early in the development of GEM-MACH. However, corrections had been made more recently to the SS emissions scheme (Sect. 3.11.3), and performance evaluation results presented in the companion paper by Moran et al. (2026) showed that the exclusion of both SS and aerosol water contributed to the consistent RAQDPS023 underprediction of PM_2.5 total mass. Sea salt was thus included in the calculation of forecast PM_2.5 total mass in RAQDPS025, although Moran et al. (2026) showed that a positive bias for SS was still present.

6.1.9 Updates to solver for gas-phase chemistry

As described in Sect. 3.4 and 3.5 both the gas-phase and aqueous-phase chemistry mechanisms in the RAQDPS023 were solved for a model chemistry time step using the computationally efficient Young and Boris (1977) predictor–corrector algorithm and employing a number of smaller integration time substeps. However, the computer code used in GEM-MACH for this solver was hardcoded, was challenging to modify and maintain, and was less accurate than some other ordinary-differential-equation (ODE) solvers now available, in particular during the day-to-night and night-to-day transitions. The publicly available Kinetic PreProcessor (KPP) software tool is a symbolic preprocessor that allows chemical mechanisms to be described in a natural and compact way (cf. Table 4), permits the user to choose from a number of different ODE solvers, and then generates computer code in three high-level computer languages (Fortran, C, Matlab) (Damian et al., 2002; Sandu and Sander, 2006). As a consequence, KPP can deliver more accurate time integrations of chemistry mechanisms, and it also reduces the effort needed to modify and maintain existing chemistry mechanisms and to implement new chemistry mechanisms.

Source code generated by the KPP to use a Rosenbrock solver was introduced in the RAQDPS025 in the gas-phase chemistry operator in place of the Young-Boris solver. One additional benefit was removal of the need to make the steady-state assumption for short-lived chemical species, which can break down in the upper troposphere due to cold temperatures.

6.1.10 Updates to representation of biogenic emissions

The representation of biogenic emissions used in the RAQDPS023 (Sect. 3.11.3) was quite dated and had an overly simple treatment of the spring and autumn transition seasons. Newer versions of the BEIS scheme were available that both provide more detailed VOC emission speciations and have updates to biomass and emission factors and to the treatment of canopy leaf temperature and soil NO emissions (https://www.epa.gov/air-emissions-modeling/biogenic-emission-inventory-system-beis, last access: 18 March 2026). Newer versions were also available for the BELD land use and vegetation database that is used in the calculation of biogenic emissions (Zhang and Moran, 2020), and many satellite-based data sets are now available that can provide high-resolution descriptions of temporal and spatial variations in vegetation phenology, including LAI (e.g., Zhang et al., 2021). For example, Campbell et al. (2022) have described the impact of changing from BEIS version 3.14 coupled with the BELD3 vegetation database to BEIS version 3.61 coupled with the BELD5 vegetation database on operational U.S. NAQFC system forecasts.

For the RAQDPS025 the BEIS scheme was updated from version 3.09 to version 3.7, the BELD vegetation database was updated from version 3 to version 4 (Zhang and Moran, 2020), and the satellite-derived, high-spatial resolution LAI monthly climatology described by Zhang et al. (2021) was also adopted, but the crude five-category treatment of phenological season (Sect. 3.9) was not changed. It is also worth noting that considerable uncertainty remains around the treatment of biogenic emissions in general (e.g., Hanna et al., 2005; Warneke et al., 2010; Guenther et al., 2012), so this component of the forecast system is likely to undergo further updates in the future.

6.1.11 Updates to anthropogenic emission inventories and emissions processing methods

The SET4.0.0 emission files used by the RAQDPS023 were based on a 2020 projected Canadian anthropogenic national emissions inventory and 2023 projected U.S. and Mexican anthropogenic national emissions inventories (Sect. 3.11.1). For the RAQDPS025 the SET4.0.0 emission files were replaced by a new set of input emission files called SET5.0.0. The SET5.0.0 emission files are based on the same U.S. and Mexican projected 2023 national inventories (just one year in time before the 2024 implementation of the RAQDPS025) as the SET4.0.0 emission files, but a 2023 projected Canadian national emissions inventory based on the 2019 APEI (ECCC, 2021) was used in place of the 2020 projected Canadian emissions inventory based on the 2015 APEI (ECCC, 2016). As was the case for the national inventories used for SET4.0.0 (Sect. S2), some pre-processing was also required for this new Canadian inventory before the emissions processing step (CMC-RAQDPS-025, 2024). Some improvements were also introduced to the emissions processing methods used by the SMOKE system to prepare the SET5.0.0 emissions files (Sect. S2). These included updates to 15 Canadian gridded spatial surrogates for several important emission sectors, including RWC (CMC-RAQDPS-025, 2024).

6.2 Future work

Work is already underway on the next version of the RAQDPS, including migration to a new upgrade of the ECCC HPC system. Possible future upgrades to GEM-MACH include improvements to the chemistry and dry deposition parameterizations, anthropogenic emissions and BB and other natural emissions, and chemical lateral boundary conditions, and implementation of chemical data assimilation.

Some improvements to the representation of atmospheric chemistry might include the addition of a parameterization of marine halogen chemistry, which can influence surface ozone concentrations in coastal regions (e.g., Sarwar et al., 2015; Li et al., 2019). Newer aqueous-phase chemistry schemes are available that account for the role of individual base cations and sea salt. Some of these schemes also include oxidation pathways for SOA formation in cloud droplets, such as organic acids, glyoxal, and isoprene epoxydiols (e.g., Gong et al., 2011; Marais et al., 2016; Fahey et al., 2017; Lamkaddam et al., 2021; Luu et al., 2025). Newer versions of ISORROPIA and other inorganic heterogeneous parameterization schemes are also available that treat individual base cations explicitly (e.g., Meng et al., 1995; Fountoukis and Nenes, 2007; Miller et al., 2024). A new solver for the aqueous-phase chemistry mechanism similar to the one described in Sect. 6.1.9 could also be implemented.

It was also noted in Sect. 3.9 and 3.11.3 that the parameterizations for gas-phase dry deposition and biogenic emissions used by the RAQDPS023 rely on a crude description of phenological season with low resolution of associated temporal and spatial variations. Many satellite-based data sets are now available that can provide high-resolution descriptions of temporal and spatial variations in vegetation phenology, including longitudinal and elevation variations (e.g., Zhang and Moran, 2020; Zhang et al., 2021). One such data set was adopted for use by the RAQDPS025 (Sect. 6.1.7 and 6.1.10), but further use of such data sets could be explored. In addition, two process terms for soil wetness not currently considered by the GEM-MACH gas-phase dry deposition scheme should be included (Sect. 3.9), and different treatments of some resistance terms could be adopted (e.g., Gao and Wesely, 1995; Wesely et al., 2002).

Regarding emissions updates, as noted in Sects. 3.11.1 and 6.1.11 there is an ongoing need to make use of the most recent anthropogenic emissions inventories to prepare updated model input emission files. Thus, the SET5.0.0 emission files used by the RAQDPS025, which have a nominal base year of 2023, could be replaced by a newer set of input emission files. Also, natural windblown dust emissions were not considered by the RAQDPS023 (Sect. 3.11.3); the addition of an aeolian dust emissions scheme would provide representation of a source of emissions that is currently missing. The estimation of BB emissions by CFFEPS might be improved by updates to fuel- and location-specific emission factors based upon Canada-specific field measurement campaigns, and recent literature (e.g., Hayden et al., 2022, Binte Shahid et al., 2024; Flood et al., 2025). The companion paper by Moran et al. (2026) also found from evaluation of RAQDPS023 forecasts that the scheme used by the RAQDPS023 to modulate fugitive dust emissions according to current meteorological conditions (Sect. 3.12) should be improved.

While the treatment of CLBCs was updated in the RAQDPS025 (Sect. 6.1.5), this updated treatment is still climatological or static in nature. In principle, CLBCs should be hour- and day-specific (e.g., Giordano et al., 2015; Ma et al., 2021; Tang et al., 2021). For example, Pendlebury et al. (2018) have described an approach whereby time-varying O₃ CLBCs can be obtained from a monthly global O₃ climatology combined with GEM model predictions of dynamic tropopause height at each time step. Another possibility to improve the CLBCs used by the RAQDPS would be to run GEM-MACH in a global configuration and then extract regional CLBCs from its forecasts (Chen et al., 2020).

Lastly, the introduction of meteorological data assimilation for the initialization of NWP models resulted in a step change in weather forecasting skill (e.g., Simmons and Hollingsworth, 2002; Bauer et al., 2015). Chemical data assimilation (CDA) offers the possibility of similar improvements for chemical weather forecasting. Some global NWP models and CTMs have implemented data assimilation of select chemical fields such as O₃, CO, NO₂, HCHO, PM_2.5, and aerosol optical depth (AOD) to improve their initialization of the chemical state (e.g., Pierce et al., 2007; Benedetti et al., 2009; Saide et al., 2013; Inness et al., 2015, 2019; Chai et al., 2017; Randles et al., 2017; Tang et al., 2017; Kumar et al., 2019; Menut and Bessagnet, 2019; Rochon et al., 2019; CMC-GDPS-8.0.0, 2021a; Ma et al., 2021). At the same time the number of available AQ measurements for North America, especially from satellite-based sensors but also from surface-based sensors, continues to increase (e.g., Naeger et al., 2021; Jaffe et al., 2023).

An important development goal for the RAQDPS is to implement an operational CDA capability, and some options to do this are available now. These include the assimilation of surface measurements based on the RDAQA system (Sect. 5.2; Robichaud, 2017; Ménard, 2021) and the use of the assimilated O₃ field from GDPS 9.0.0 forecasts (Sect. 6.1.2; Rochon et al., 2019; CMC-GDPS-8.0.0, 2021a).

7 Summary

Environment and Climate Change Canada (ECCC) has run an operational, continental-scale Regional Air Quality Deterministic Prediction System (RAQDPS) for North America since 2001. The earliest version of the RAQDPS employed an offline chemical transport model called CHRONOS that used the ECCC Global Environmental Multiscale (GEM) numerical weather prediction model as a meteorological driver. In 2009 CHRONOS was replaced in the RAQDPS by the online, one-way-coupled GEM-MACH (Global Environmental Multiscale–Modeling Air quality and CHemistry) chemical weather model, which consists of GEM as a meteorological host model and an embedded chemistry module named MACH (Modelling Air quality and CHemistry) that can be activated (or not) by a namelist setting. Since 2009 there have been 24 upgrades of varying magnitude made to the GEM-MACH-based version of the RAQDPS to improve forecast skill and to adjust to changes in upstream meteorological systems and ECCC's supercomputers (Table A1). This paper describes version 023 of the RAQDPS, including its dynamics, physics, chemistry, numerics, inputs, outputs, related upstream meteorological systems and downstream AQ systems, and recent and potential future developments. Significantly, it is the first comprehensive publication describing any RAQDPS version. A companion paper (Moran et al., 2026) provides a comprehensive evaluation of RAQDPS023 performance for its first forecast year and four hindcast years.

The RAQDPS023 (or RAQDPS-OP023) became the ECCC operational regional AQ forecast system on 1 December 2021; it was then migrated without any algorithmic changes to a new ECCC computer system on 28 June 2022 and renamed RAQDPS-OP024. It then continued to run operationally for another two years until it was replaced by an upgraded version, the RAQDPS025, on 11 June 2024. A special wildfire version of the RAQDPS-OP023 called FireWork023 or RAQDPS-FW023, which was identical to the RAQDPS-OP023 except for the addition of near-real-time (NRT) biomass burning (BB) emissions as an input, went operational at the same time, was migrated to the new computer system at the same time (RAQDPS-FW024), and was replaced at the same time in 2024. Running this second member of the ECCC regional AQ forecast suite allowed forecasts of the transport and diffusion of wildfire smoke plumes to be calculated in a post-processing step as the difference between the PM_2.5 forecasts from the two model versions.

The RAQDPS-OP023 and RAQDPS-FW023 both used GEM-MACH version 3.1.0.0, which was built on GEM version 5.1.0. Both systems employed a regional configuration with 10 km horizontal grid spacing and a geographic domain covering Canada, the continental U.S. and most of Alaska, and northern Mexico. Both systems were run twice each day starting at 00:00 and 12:00 UTC to produce 72 h forecasts. The RAQDPS-FW023 obtained NRT satellite-detected BB emission forecasts for Canada and the U.S. for the next 72 h from version 4.1 of the Canadian Forest Fire Emissions Prediction System (CFFEPS) before each run (Sect. S3).

To obtain the meteorological initial conditions and boundary conditions that were needed, each RAQDPS023 run followed a run for the same forecast period by the corresponding ECCC operational regional weather forecast system, the 10 km Regional Deterministic Prediction System 8.0.0. The RDPS 8.0.0 was a regional configuration of GEM 5.1.0 that made 84 h forecasts four times per day. Its horizontal grid was a superset of the RAQDPS023 grid. To obtain the meteorological initial conditions and boundary conditions that it needed in turn, each RDPS 8.0.0 run followed a run for the same forecast period by the ECCC operational global weather forecast system, the 15 km Global Deterministic Prediction System 8.0.0. The GDPS 8.0.0 was a global configuration of GEM 5.1.0 that made 144 h forecasts four times per day. For maximum consistency the GDPS 8.0.0, RDPS 8.0.0, and RAQDPS023 all used the same rotated latitude-longitude map projection and the same 84 staggered hybrid vertical levels with a top at 0.1 hPa.

All three systems also employed the same GEM dynamical core, including the hydrostatic assumption. Time integration was performed using a two-time-level, iterative-implicit, three-dimensional semi-Lagrangian scheme. The GDPS 8.0.0 used a 450 s dynamics integration time step while both the RDPS 8.0.0 and RAQDPS023 used a 300 s time step. Since semi-Lagrangian advection schemes do not conserve mass, a shape-preserving monotonicity correction and a global mass fixer were both used for tracer advection for all three systems. Mass conservation is very relevant for a regional AQ model since there can be many sharp chemical concentration gradients in the vicinity of strong local emissions sources, and the resulting loss of shape preservation after advection can introduce physically unrealistic negative concentrations or positive overshoots. The same physics parameterizations were also used by the three systems (Table 1), with only very minor differences in a few parameter settings between the global and regional versions to account for their different horizontal grid spacing. And since GEM-MACH is an online model the RAQDPS023 could access many meteorological variables predicted by GEM 5.1.0 that were needed by some MACH chemistry parameterizations (Table 2).

The RAQDPS023 predicted 41 gas-phase species and 18 particle-phase chemical species, but its three key AQ predictands were NO₂, O₃, and PM_2.5 total mass, which are needed to calculate the Canadian Air Quality Health Index. The AQHI is a health-based, multi-pollutant, additive, no-threshold hourly AQ index that is used to communicate current and predicted AQ levels to the Canadian public. Of these three predictands, PM_2.5 mass poses a particular challenge because it is a complex pollutant: it has a wide size range; it is composed of many elements and compounds, both inorganic and organic; it is subject to many physical and chemical processes; and it has many types of emissions sources.

The MACH chemistry module employs a sectional representation to model the PM size distribution, where the number of size bins can be specified by the user. For the RAQDPS023 as well as earlier versions, a major simplification made to ensure that system execution times would be short enough to fit within its operational 30 min time slot was to consider just two PM size bins, one corresponding to particle Stokes diameters in the 0–2.5 µm range (“fine fraction” or PM_2.5) and the other to the 2.5–10 µm range (“coarse fraction” or PM_cf). Together these two size bins constitute PM₁₀. This is a very coarse representation of the PM size distribution, but it still allowed the model to predict concentrations of PM_2.5 and PM₁₀, which are the two PM size ranges considered by most AQ measurement networks and emission inventories. MACH for the RAQDPS023 also represented PM_2.5 chemical composition with a limited number of chemical components, namely sulfate, nitrate, ammonium, elemental carbon, primary organic matter, secondary organic matter, crustal material, sea salt, and aerosol water. However, these nine chemical components align fairly well with the chemical components that are measured by PM_2.5 speciation networks, thus permitting evaluation of PM_2.5 composition predictions. A third simplification in the treatment of PM was that each size bin was assumed to be internally mixed, an assumption that is the least good close to sources of PM but which improves as particles age in the atmosphere.

MACH contains parameterizations to represent a number of physical and chemical processes that influence PM mass concentrations. As described in Sect. 3 these include anthropogenic emissions, smokestack plume rise, sea-salt and biogenic emissions, particle nucleation, condensational growth, coagulation, gravitational settling and dry deposition, hygroscopic growth and aerosol particle activation, vertical diffusion, gas-phase chemistry, secondary organic aerosol formation, gas-particle partitioning, inorganic heterogeneous chemistry, aqueous-phase chemistry, and in-cloud and below-cloud precipitation scavenging.

Chemical tracer initial conditions for the RAQDPS023 were copies of the chemical tracer fields corresponding to that initial time from the previous forecast, a so-called perpetual forecast approach. Chemical lateral boundary conditions were obtained from a set of “climatological” seasonal vertical cross-sections of model species mass mixing ratios (MMRs) supplied for each lateral boundary. One other simplification made for the RAQDPS023 to reduce execution time was to call the MACH module every third dynamics time step, that is, to employ a chemistry time step of 900 s. This simplification was possible due to the assumption that there was no feedback of chemistry on meteorology. Even so, activation of the MACH chemistry module increased RAQDPS023 execution time by a factor of 4.4 compared to running only the GEM code.

The pre-processed, hourly, gridded anthropogenic emissions input files used by the RAQDPS-OP023 and RAQDPS-FW023, which are referred to as SET4.0.0 emissions, were new to this version (Table A1). They were based on three national emission inventories: (i) a projected 2020 Canadian Air Pollutant Emissions Inventory; (ii) a projected 2023 U.S. National Emissions Inventory (NEI); and (iii) a projected 2023 Mexican NEI. The SET4.0.0 files were generated during the development of the RAQDPS-OP023 and well before its operational deployment using version 4.7 of the SMOKE emissions processing system. Hourly BB emissions, on the other hand, were forecast for the next 72 h in advance of each RAQDPS-FW023 run by CFFEPS v4.1. Biogenic emissions were treated in a third way; they were estimated every chemistry time step by inline code in MACH. Sea-salt emissions were also estimated every chemistry time step.

The RAQDPS023 routinely output a large number of hourly gridded fields, including 41 gas-phase chemical MMR fields, eight PM_2.5 and eight PM_cf chemical component MMR fields, PM_2.5 and PM₁₀ total mass concentration fields, 21 dry deposition flux fields, 12 wet deposition flux fields, and four biogenic emission flux fields. Two downstream operational systems made use of some of these outputs. Point-specific hourly forecasts of NO₂, O₃, and PM_2.5 concentrations at Canadian AQ measurement station locations were fed to the Updateable Model Output Statistics for Air Quality (UMOS-AQ) system. UMOS-AQ is a statistical post-processing package for bias correction that can compensate for systematic model errors and account for unresolved SGS phenomena. Its hourly, station-specific regression coefficients are recalculated weekly. A second post-processing system, version 2.0.0 of the Regional Deterministic Air Quality Analysis (RDAQA) system, used RAQDPS023 hourly gridded forecast fields as first-guess fields together with surface AQ measurements to generate objective analyses of hourly North American surface concentration fields for O₃, NO₂, NO, PM_2.5, PM₁₀, SO₂, and the AQHI. RAQDPS023 forecast performance was also monitored on an ongoing basis, including by a quarterly comparison with AQ forecasts for North America made by four international operational models.

The paper concludes by describing the main innovations made for the RAQDPS025, the successor to the RAQDPS023, which was implemented operationally in June 2024. One key innovation of the RAQDPS025 was to combine the RAQDPS-OP and RAQDPS-FW systems into a single unified system while still identifying wildfire smoke plumes. The RAQDPS025 also has updates to the treatments of vertical diffusion, gas-phase chemistry, gas-phase dry deposition, below-cloud scavenging of size-resolved PM, chemical lateral boundary conditions, and biogenic emissions. As well the meteorological piloting model was changed (to GDPS 9.0.0) and a new set of anthropogenic input emissions files (SET5.0.0) was introduced. Lastly, a short summary of some possible areas of improvements to future RAQDPS versions is provided.

Appendix A: List of acronyms and abbreviations

ADAGIO	Atmospheric Deposition Analysis Generated by optimal Interpolation from Observations
ADOM	Acid Deposition and Oxidant Model
AGL	above ground level
ALKA	See Table 3
ALKE	See Table 3
AM	ammonium
AOD	aerosol optical depth
APEI	Air Pollutant Emissions Inventory (Canada)
AQ	air quality
AQHI	Air Quality Health Index (Canada)
ARCTAS	Arctic Research of the Composition of the Troposphere from Aircraft and Satellites
AROM	See Table 3
AURAMS	A Unified Regional Air quality Modelling System
BB	biomass burning
BC	black carbon
BEIS	Biogenic Emission Inventory System
BELD3	Biogenic Emissions Landuse Database version 3
CAC	criteria air contaminant
CAM	Canadian Aerosol Module
CAM-chem	Community Atmosphere Model with chemistry
CAPE	convective available potential energy
CCN	cloud condensation nuclei
CDA	chemical data assimilation
CF	coarse fraction
CFC	chlorofluorocarbon
CFFEPS	Canadian Forest Fire Emissions Prediction System
CFS	Canadian Forest Service
CHRONOS	Canadian Hemispheric and Regional Ozone and NOx System
CLBC	chemical LBC
CM	crustal material
CMAQ	Community Multiscale Air Quality model

CMAS	Community Modeling and Analysis System (University of North Carolina)
CMC	Canadian Meteorological Centre
CS	condensable species
CSN	Chemical Speciation Network (U.S.)
CTM	chemical transport model
CWFIS	Canadian Wildland Fire Information System
DDF	dry deposition scaling factor
DOW	day of the week
EC	elemental carbon
ECCC	Environment and Climate Change Canada
EOTH	emissions of “other” VOCs (NAPAP classes 21, 23, 25, 26, 31, and 32)
EPA	Environmental Protection Agency (U.S.)
FBP	(Canadian) Forest Fire Behaviour Prediction subsystem
FEPS	Fire Emission Production Simulator (U.S.)
FW	FireWork
FWI	(Canadian) Forest Fire Weather Index subsystem
GAFIS	Global Air quality Forecasting and Information System initiative (WMO)
GDPS	Global Deterministic Prediction System
GEM	Global Environmental Multiscale (model)
GEM-MACH	Global Environmental Multiscale–Modelling Air quality and CHemistry (model)
GMM	gram molar mass
HETV	HETerogeneous Vectorized scheme
HPC	high-performance computing
IAU	incremental analysis update
IAY	instantaneous aerosol yield
IC	initial condition
ILMC	iterative locally mass conserving (scheme)
IMPROVE	Interagency Monitoring of Protected Visual Environments network (U.S.)
ISBA	Interactions between Soil–Biosphere–Atmosphere (land surface model)
ISOP	isoprene (see Table 3)
ISORROPIA	“equilibrium” (modern Greek)
KPP	Kinetic PreProcessor
LAI	leaf area index

LAM	limited area model
LBC	lateral boundary condition
LINOZ	LINearized OZone mechanism
LST	local standard time
LWC	liquid water content
MACH	Modelling Air quality and CHemistry (module)
MLR	multi-variate linear regression
MMR	mass mixing ratio
MODIS	Moderate Resolution Imaging Spectroradiometer
MOS	model output statistics
MOZART4	Model for OZone And Related chemical Tracers version 4
MPI	Message Passing Interface
NAPAP	National Acid Precipitation Assessment Program (U.S.)
NAPS	National Air Pollution Surveillance system (Canada)
NAQFC	National Air Quality Forecast Capability (U.S.)
NASA	National Aeronautics and Space Administration (U.S.)
NEI	National Emissions Inventory (U.S., Mexico)
NI	nitrate
NMHC	non-methane hydrocarbon
NOAA	National Oceanic and Atmospheric Administration (U.S.)
NRT	near-real time
NWP	numerical weather prediction
OA	objective analysis
ODE	ordinary differential equation
OH	hydroxl radical
OI	optimal interpolation
OpenMP	Open Multi-Processing
OVOC	other VOCs
PAR	photosynthetic active radiation
PBL	planetary boundary layer
PM	particulate matter
PMcf	PM coarse fraction (2.5–10 µm)
PM2.5	PM with aerodynamic diameter of less than 2.5 µm
PM10	PM with aerodynamic diameter of less than 10 µm
POM	primary organic matter
QA/QC	quality assurance/quality control
RAQDPS	Regional Air Quality Deterministic Prediction System
RAQDPS-FW	Operational RAQDPS with FireWork BB emissions
RAQDPS-OP	Operational RAQDPS without FireWork BB emissions
RDAQA	Regional Deterministic Air Quality Analysis

RDPS	Regional Deterministic Prediction System
RH	relative humidity
ROS	rate of spread
RVFSP-WAC	North American Regional Vegetation Fire and Smoke Pollution Warning and Advisory Centre
RWC	residential wood combustion
SCC	source classification code
SGS	subgrid-scale
SI	Système International d'unités
SMOKE	Sparse Matrix Operation Kernel Estimation system (CMAS)
SOA	secondary organic aerosol
SOM	secondary organic matter
SS	sea salt
SU	sulfate
TF	transportable fraction
TKE	turbulent kinetic energy
TOLU	See Table 3
UMOS-AQ	Updateable Model Output Statistics–Air Quality
USFS	U.S. Forest Service
UTC	Universal Time Coordinated
VIIRS	Visible Infrared Imaging Radiometer Suite
VIT	vehicle-induced turbulence
VOC	volatile organic compound
WA	aerosol water
WHO	World Health Agency
WMO	World Meteorological Organization
4D-IAU	four-dimensional incremental analysis update

List of symbols

d	particle diameter (cm)
D	raindrop or frozen hydrometeor diameter (cm)
Di	gas-phase diffusivity of species i (cm2 s−1)
Dm	mean raindrop or frozen hydrometeor diameter (cm)
$\overline{E}(d,l)$	mean snowflake-particle collection efficiency (unitless)
E(d,Dm)	raindrop-particle collection efficiency (unitless)
$\overline{E}(d,D_{\mathrm{m}})$	mean raindrop-particle collection efficiency (unitless)
fcld	total cloud cover fraction by layer (unitless)
fctr	precipitation production rate (s−1)
H′	effective equilibrium constant for gas-liquid partitioning
K	equilibrium constant for aqueous-phase reactions
Kom,1i	first partitioning coefficient of the ith condensable species (µg m−3)
Kom,2i	second partitioning coefficient of the ith condensable species (µg m−3)
Kz	vertical eddy diffusivity for tracers (m2 s−1)
L	characteristic collection length scale of frozen hydrometeor (may not be D) (cm)
L	liter
Lw	volumetric liquid water fraction in air (m3 H2O per m3 air)
N	number of size bins or sections
Nd	cloud droplet number concentration (cm−3)

Np	particle number density (cm−3)
Np,i	particle number density in size bin i (cm−3)
P	(liquid) precipitation rate (mm h−1)
P	atmospheric pressure (hPa)
ps	surface pressure (hPa)
R	universal gas constant (J K−1 mol−1 or L atm K−1 mol−1)
Ra	aerodynamic resistance (s m−1)
Rb	quasi-laminar sublayer resistance (s m−1)
Rc	surface resistance (s m−1)
Ri	gradient Richardson number (unitless)
Rib	bulk Richardson number (unitless)
T	dry-bulb temperature (K)
u∗	friction velocity (m s−1)
Vrt	raindrop terminal fall speed (m s−1)
Vst	frozen hydrometeor terminal fall speed (m s−1)
α1i	mass-based stoichiometric coefficient 1 for SOA formation from ith CS (unitless)
α2i	mass-based stoichiometric coefficient 2 for SOA formation from ith CS (unitless)
λ	turbulence-regime-specific mixing length (m)
λe	dissipation length scale (m)
Λr(d)	rain scavenging rate coefficient for aerosol particles of diameter d (s−1)
Λs(d)	snow scavenging rate coefficient for aerosol particles of diameter d (s−1)
Λi	rain scavenging rate coefficient for ith gas (s−1)

Table A1 Chronology of operational RAQDPS and RAQDPS-FW versions. Boldface font indicates a major update to the RAQDPS version and the RAQDPS-FW version (which was introduced operationally in April 2016 as a seasonal forecast system). (See also ECCC changelogs for the RAQDPS (https://eccc-msc.github.io/open-data/msc-data/nwp_raqdps/changelog_raqdps_en/, last access: 18 March 2026) and RAQDPS-FW (https://eccc-msc.github.io/open-data/msc-data/nwp_raqdps-fw/changelog_raqdps-fw_en/, last access: 18 March 2026).)