Description of the uEMEP_v5 downscaling approach for the EMEP MSC-W chemistry transport model

A description of the new air quality downscaling model – the urban EMEP (uEMEP) and its combination with the EMEP MSC-W model (European Monitoring and Evaluation Programme Meteorological Synthesising Centre West) – is presented. uEMEP is based on well-known Gaussian modelling principles. The uniqueness of the system is in its combination with the EMEP MSC-W model and the “local fraction” calculation contained within it. This allows the uEMEP model to be imbedded in the EMEP MSC-W model and downscaling can be carried out anywhere within the EMEP model domain, without any double counting of emissions, if appropriate proxy data are available that describe the spatial distribution of the emissions. This makes the model suitable for high-resolution calculations, down to 50 m, over entire countries. An example application, the Norwegian air quality forecasting and assessment system, is described where the entire country is modelled at a resolution of between 250 and 50 m. The model is validated against all available monitoring data, including traffic sites, in Norway. The results of the validation show good results for NO₂, which has the best known emissions, and moderately good for PM₁₀ and PM_2.5. In Norway, the largest contributor to PM, even in cities, is long-range transport followed by road dust and domestic heating emissions. These contributors to PM are more difficult to quantify than NO_x exhaust emission from traffic, which is the major contributor to NO₂ concentrations. In addition to the validation results, a number of verification and sensitivity results are summarised. One verification showed that single annual mean calculations with a rotationally symmetric dispersion kernel give very similar results to the average of an entire year of hourly calculations, reducing the runtime for annual means by 4 orders of magnitude. The uEMEP model, in combination with EMEP MSC-W model, provides a new tool for assessing local-scale concentrations and exposure over large regions in a consistent and homogenous way and is suitable for large-scale policy applications.

1 Introduction

The EMEP MSC-W model is a chemistry transport model which has been developed by the Meteorological Synthesizing Centre – West (MSC-W) of EMEP, the European Monitoring and Evaluation Programme, under the UN Convention on Long-range Transboundary Air pollution (LRTAP). It is documented in Simpson et al. (2012) and has been used for many years mainly for regional but also for global applications. The aim of uEMEP (urban EMEP) is to further extend the application of the EMEP MSC-W chemical transport model down to near-street-level modelling. The downscaling methodology builds on classical Gaussian plume modelling and is integrated with the EMEP MSC-W model physical parameterisations and emission data in such a way as to provide a consistent model description from regional to local scales. Unlike other urban-scale models, uEMEP is intended not just to achieve local-scale modelling for an individual city or area but to provide local-scale modelling over entire countries or continents, providing high-resolution modelling over large areas and allowing air quality assessment and exposure calculations at high resolution everywhere. Air quality modelling is often separated into global, regional, urban and local scales. In this context, “local” refers to the ability of the model not just to represent urban background concentrations but also to represent street-level concentrations. There are a number of models that span the global or regional scale where grid resolutions down to 4–10 km have been achieved, e.g. EMEP MSC-W (Simpson et al., 2012; Werner et al., 2018), CHIMERE (Menut et al., 2013), SILAM (Sofiev et al., 2015), LOTOS-EUROS (Kranenburg et al., 2013), MATCH (Andersson et al., 2007) and CMAQ (Appel et al., 2017). There are a number of Gaussian modelling systems that cover the urban and local scales over limited areas, usually individual cities, e.g. ADMS (Stocker et al., 2012) and AERMOD (Cimorelli et al., 2004). In addition, there are some limited area models that combine Eulerian and Gaussian plume type models in a single system, e.g. Karamchandani et al. (2009), Kim et al. (2018) and Karl et al. (2019). If the full model cascade is to be achieved, such as the THOR forecast system in Denmark (Brandt et al., 2001), then this requires linking different model systems together to achieve high-resolution calculations in a limited area. An alternative approach to achieving high-resolution concentration fields over large regions without the use of chemistry transport models (CTMs) (Chemical transport models) are land use regression methods (e.g. Vizcaino and Lavalle, 2018); however, their lack of underlying physics does not make them useful for planning or policy applications.

Earlier work on the downscaling of CTMs over large regions include the population covariance correction factor from Denby et al. (2011) and the dispersion kernel methods from Theobald et al. (2016) and Maiheu et al. (2017). There are similarities between uEMEP and these last two studies as both use Gaussian models for the downscaling. The major difference between uEMEP and these two Gaussian kernel methods is that uEMEP can be applied on hourly data, as well as annual data, and that uEMEP is integrated with the “local fraction” scheme in EMEP MSC-W (Wind et al., 2020) that avoids double counting of emissions in a consistent manner.

In this paper, we provide an overall description of the uEMEP methodology and how it is combined with the “local fraction” scheme in EMEP MSC-W (Sect. 2). The uEMEP model physical parameterisations are then given in Sect. 3. In Sect. 4, an application example of the methodology, the Norwegian air quality forecasting service, is described. Validation of the modelling system against all Norwegian monitoring data is presented in Sect. 5 together with a summary of verification and sensitivity studies. Various aspects of the modelling are discussed in Sect. 6 and conclusions made in Sect. 7. More detailed information on the parameterisations, validation and verification is also provided in the Supplement.

2 Methodology

In this section, we describe the concepts and methodologies for the application of the coupled modelling system uEMEP and EMEP MSC-W.

2.1 Overall concept and methodology

uEMEP provides a consistent methodology for redistributing/downscaling gridded CTM emissions and concentrations, in this case from the EMEP MSC-W model, to higher-resolution subgrids within the CTM grids. Proxy data, that represent the spatial distribution of the emissions, are used to redistribute emissions to subgrids. Concentrations are then recalculated at the subgrid level, using a Gaussian model, whilst removing the local contribution of the CTM at these subgrids but keeping the non-local contributions. This procedure consistently avoids double counting of emissions.

Throughout this paper, we refer to the downscaling grids in uEMEP as “subgrids”. These may be any size but current applications use a range of between 25 and 250 m. When referring to the EMEP MSC-W model, we use the term “grid”. This may also vary dependent on the application but is usually in the range of 2 to 15 km. The term “local” is also used. Local, in regard to EMEP, means the local EMEP grids, so terms such as “local fraction” refer to a particular grid and the other EMEP grids in the “local region”, for example, within a range of 5×5 EMEP grids. When referring to “local” in uEMEP, we also refer to subgrid contributions from within this local EMEP region. This is visualised in Fig. 1. “Non-local”, in regard to uEMEP, refers to any contribution that is not downscaled or calculated with uEMEP, usually contributions from outside the local EMEP region but these can also be other source sectors not accounted for by uEMEP. We will refer to concentrations provided by the EMEP MSC-W model simply as EMEP concentrations.

uEMEP can be run using two downscaling methods, both of which make use of Gaussian dispersion modelling to describe the redistribution of concentrations at high resolution. Both methods make use of the recently developed “local fraction” functionality in the EMEP model (Sect. 2.3; Wind et al., 2020) to avoid double counting of emissions and to allow near-seamless integration of the two models. The two downscaling methods are as follows:

1. Emission redistribution. Redistribution of the EMEP gridded emissions using emission proxy data to subgrids and calculation of concentrations through Gaussian modelling, with removal of the locally emitted EMEP concentrations.

2. Independent emissions. Independent high-resolution emission data on subgrids and calculation of the concentrations through Gaussian modelling, with removal of the locally emitted EMEP source contributions. The independent emissions do not need to be consistent with the EMEP gridded emissions in this case.

In addition, calculations can be made on either long-term mean emissions, using a rotationally symmetric dispersion kernel (Sect. 3.2), or on hourly emissions, using a slender plume Gaussian dispersion model (Sect. 3.1).

Typical source sectors downscaled using uEMEP include traffic, residential combustion, shipping and industry. The sectors addressed will depend on the availability of high-resolution data for distributing them. uEMEP is only applied to the primary emissions of PM₁₀, PM_2.5 and NO_x and does not involve any complex chemistry or secondary formation of particles. The concentrations of NO₂ and O₃ are calculated with uEMEP using a simplified chemistry scheme (Sect. 3.4 and 3.5).

2.2 Subgrid calculation method

The choice of downscaling method will depend on the quality of the high-resolution proxy or emission data available, whether the calculation will be made on hourly or annual data and on the EMEP model resolution. The first downscaling method, emission redistribution, will be applied when only approximate proxy data for emissions are available, which will be the case in many large-scale applications. Examples of such proxy data may be population density, road network data or land use data. The second downscaling method, independent emissions, is suitable for when high-quality emission data are available that are consistent between uEMEP and EMEP. When the gridded emission data are entirely consistent with the proxy data, i.e. the proxy data are given as emissions and are aggregated to the CTM grid emissions, the two methods are equivalent.

The total subgrid concentrations C_SG(i,j) at subgrid indexes (i,j) are calculated by adding the local Gaussian concentration $C_{\mathrm{SG},\mathrm{local}}(i,j)$ and the non-local part of the EMEP grid concentration $C_{\mathrm{SG},\mathrm{nonlocal}}(i,j)$ and is written simply as

$$\begin{array}{}\text{(1)}& C_{\mathrm{SG}}\left(i,j\right)=C_{\mathrm{SG},\mathrm{local}}\left(i,j\right)+C_{\mathrm{SG},\mathrm{nonlocal}}\left(i,j\right),\end{array}$$

where we use the subscript notation “SG” to denote any subgrid value and in further equations the subscript “G” to indicate any EMEP grid value. $C_{\mathrm{SG},\mathrm{local}}(i,j)$ is determined directly from the subgrid dispersion calculation:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& C_{\mathrm{SG},\mathrm{local}}\left(i,j\right)=\sum _{i^{\prime }=\mathrm{1}}^{n_x}\sum _{j^{\prime }=\mathrm{1}}^{n_y}{\displaystyle \frac{E_{\mathrm{SG},\mathrm{local}}\left(i^{\prime },j^{\prime }\right)}{U\left(i,j,i^{\prime },j^{\prime }\right)}}\\ & \cdot I\left(r\left(i,j,i^{\prime },j^{\prime }\right),h(i^{\prime },j^{\prime }),z\left(i,j\right)\right).\end{array}\end{array}$$

Here, E_SG,local is the emission attributed to each subgrid and U is the wind speed, which in uEMEP is dependent on both the source and receptor subgrid values (Sect. 3.1). n_x and n_y represent the extent of the subgrid calculation window. The function $I(r,h,z)$ is the dispersion intensity function (Sect. 3.1) that determines the dispersion of the emission source E_SG,local with the horizontal spatial vector $r(i,j,i^{\prime },j^{\prime })$ between the receptor grid points (i,j) and the source grid points ($i^{\prime },j^{\prime }$) at height z(i,j). The source height $h(i^{\prime },j^{\prime })$ is also specified. The contribution from every proxy emission subgrid ($i^{\prime },j^{\prime }$), within the defined subgrid calculation window (n_x,n_y) is calculated and summed at the receptor subgrid (i,j) centred within subgrid calculation window; see Fig. 1.

When using the emission redistribution method, E_SG,local is calculated using the EMEP grid emissions E_G(I,J) and the proxy data for emissions, P_emission. P_emission is normalised within the EMEP grid in the following way to determine the subgrid emission E_SG,local:

$$\begin{array}{}\text{(3)}& E_{\mathrm{SG},\mathrm{local}}\left(i^{\prime },j^{\prime }\right)=E_{\mathrm{G}}(I,J){\displaystyle \frac{P_{\mathrm{emission}}\left(i^{\prime },j^{\prime }\right)}{\sum _{i^{\prime }=\mathrm{1}}^{n_x}\sum _{j^{\prime }=\mathrm{1}}^{n_y}{P}_{\mathrm{emission}}\left(i^{\prime },j^{\prime }\right)}}.\end{array}$$

When using the independent emission method, the local subgrid emissions E_SG,local are specified directly.

2.3 Calculation of the non-local contribution from EMEP

The term $C_{\mathrm{SG},\mathrm{nonlocal}}(i,j)$, given in Eq. (1) and further derived in Eqs. (9) and (10), is the non-local contribution from the EMEP grid at the specific subgrid (i,j). Though this is based on the non-local contribution provided by EMEP at grids (I,J) interpolation due to the moving window (Sect. 2.4) surrounding each receptor subgrid means that non-local contributions are specified at the subgrid level. The gridded non-local contribution, $C_{\mathrm{G},\mathrm{nonlocal}}(I,J)$, is derived from the “local fraction” calculation in EMEP. The methodology is described completely in Wind et al. (2020) but the essential elements are reproduced here.

The local fraction methodology corresponds to a tagging method where pollutants from different origins are tagged and stored individually. In this case, the tagging occurs relative to the surrounding grid cells of any individual grid. This means that emissions from any grid cell are tagged and followed through the various model processes out to neighbouring grid cells. Tagged species are assumed to be inert species, primary PM and NO_x, for this downscaling application as chemical reactions are not included in the tagging. It is generally not computationally possible, or in this application necessary, to follow all grid cell contributions to all other grids within the EMEP model domain. The local fraction region extent is then limited. In Wind et al. (2020), the local fraction region extent (n_lf) was tested up to a 161×161 EMEP grids on low-resolution EMEP runs for Europe. Generally, 21×21 EMEP grids were found to be computationally and memory efficient. In the uEMEP application, the local fraction region needs only be as large as the uEMEP calculation window, i.e. the allowed distance from the receptor subgrid to the emission subgrids. In the forecasting application discussed in Sect. 4, this requires only a 5×5 EMEP grid local fraction region. Sensitivity to the size of this region is discussed in Sect. S5.2. For each EMEP grid (I,J), there will be an associated local fraction grid LF($I,J,I_{\mathrm{lf}},J_{\mathrm{lf}}$) that specifies the fraction of the surrounding grids contributing to the (I,J) grid. I_lf and J_lf are indexed from $-n_{\mathrm{lf}}/\mathrm{2}$ to $+n_{\mathrm{lf}}/\mathrm{2}$.

With use of the local fraction, the local (C_G,local) and non-local (C_G,nonlocal) contributions from any particular primary pollutant to an EMEP grid (I,J) are given by the sum of all the contributing local fraction contributions of the local sources (s=1 to n_source) and the non-local contribution, specified by

$$\begin{array}{}\text{(4)}& {\displaystyle }{C}_{\mathrm{G},\mathrm{local}}\left(I,J,I_{\mathrm{lf}},J_{\mathrm{lf}},s\right)=\mathrm{LF}\left(I,J,I_{\mathrm{lf}},J_{\mathrm{lf}},s\right)C_{\mathrm{G}}\left(I,J\right)\text{(5)}& {\displaystyle }\begin{array}{rl}& C_{\mathrm{G},\mathrm{nonlocal}}\left(I,J\right)=C_{\mathrm{G}}\left(I,J\right)\\ & -\left\{\sum _{I^{\prime }=-n_{\mathrm{mw}}/\mathrm{2}}^{+n_{\mathrm{mw}}/\mathrm{2}}\sum _{J^{\prime }=-n_{\mathrm{mw}}/\mathrm{2}}^{+n_{\mathrm{mw}}/\mathrm{2}}\sum _{s=\mathrm{1}}^{n_{\mathrm{source}}}{C}_{\mathrm{G},\mathrm{local}}\left(I,J,I^{\prime },J^{\prime },s\right)\right\}.\end{array}\end{array}$$

Note that in Wind et al. (2020) C_G,local and C_G are termed LP (local pollutant) and TP (total pollutant), respectively, and the local fraction index is specified here using (I_lf,J_lf) instead of (Δx_s⋅Δy_s). This change is for compatibility with the notation used for the uEMEP application.

2.4 Moving window calculation of local and non-local EMEP contributions

When determining the local and non-local EMEP contribution at any uEMEP subgrid receptor, a moving window methodology is applied. The aim of the moving window calculation is to represent as well as possible the local and non-local EMEP contributions at any one subgrid, in effect creating an EMEP grid that is centred on the receptor subgrid. The moving window is centred on the receptor subgrid (i,j) and its size is specified by the number of EMEP grids it covers (n_mw,n_mw). The moving window region is the same as the uEMEP calculation window in extent, which is also defined by the number of subgrids (n_x,n_y) (Sect. 2.2). n_mw is given by the user but it must not be larger than the area covered by the EMEP local fraction region (n_lf), i.e. $n_{\mathrm{mw}}\le n_{\mathrm{lf}}-\mathrm{1}$. Figure 1 shows an example where n_lf=5 and n_mw=4.

Since we need to account for all source contributions from EMEP within the moving window and since the subgrids are not centred in the middle of the EMEP grids, the local contribution from the EMEP grids for any particular source sector s can be written as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& C_{\mathrm{G},\mathrm{local}}\left(i,j,s\right)=\sum _{I^{\prime }=-n_{\mathrm{mw}}/\mathrm{2}}^{+n_{\mathrm{mw}}/\mathrm{2}}\sum _{J^{\prime }=-n_{\mathrm{mw}}/\mathrm{2}}^{+n_{\mathrm{mw}}/\mathrm{2}}\\ & C_{\mathrm{G},\mathrm{local}}\left(I,J,I^{\prime },J^{\prime },s\right)\cdot w\left(i,j,I^{\prime },J^{\prime },s\right).\end{array}\end{array}$$

Here, the weighting variable $w(i,j,I^{\prime },J^{\prime },s)$ refers to the weighting of the EMEP grid relative to the receptor subgrid and the index I,J refers to the EMEP grid which contains the uEMEP subgrid (i,j). For EMEP grids entirely within the moving window, this weighting will be unity, but for EMEP grids only partially within the moving window this weighting will be less than unity as part of that EMEP grid will also contribute to the non-local concentrations.

There are two methods implemented in uEMEP for specifying these weights. The simplest and most often used is area weighting where only the area fraction of the EMEP grid that is within the moving window for that particular receptor subgrid is included in the local contribution. This is illustrated in Fig. 1 and is usually sufficient for the calculation, especially when the number of EMEP grids covered by the moving window is larger than 3×3. Mathematically, the area weighting, w_a, can be written as

$$\begin{array}{}\text{(7)}& w_{\mathrm{a}}\left(i,j,I^{\prime },J^{\prime }\right)={\displaystyle \frac{\left\{a(i,j)\cap A(I^{\prime },J^{\prime })\right\}}{A(I^{\prime },J^{\prime })}},\end{array}$$

where $A(I^{\prime },J^{\prime })$ is the area and position of each EMEP grid, a(i,j) is the area and position of the moving window centred at the receptor subgrid point (i,j) and $a\left(ij\right)\cap A\left(I^{\prime }{J}^{\prime }\right)$ is the overlapping area of these two regions. For the case where n_mw=1, this area weighting is equivalent to a bilinear interpolation of the surround EMEP grids. Area weighting is not dependent on the source.

Figure 1 Schematic representation of the moving window region. It shows the regions used for uEMEP calculations and the area and emission weighting selection used to determine the local and non-local EMEP contributions at the calculation (receptor) subgrid. The extent of the subgrids is only partially shown.

When the moving window only covers a limited number of EMEP grids and when high-resolution emission data are used that are compatible with the EMEP grid emissions, this weighting can also be based on the high-resolution emission data. This better represents the moving window concept because it reflects the effect of moving the EMEP grid to be centred on the receptor subgrid in a more realistic way. In this case, the emission weighting term (w_e) on the edge of the moving window will be determined by the fraction of the total subgrid emissions within the moving window and within the EMEP grid, instead of the area. This can be written as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& w_{\mathrm{e}}\left(i,j,I^{\prime },J^{\prime },s\right)=\\ & {\displaystyle \frac{\sum e\left(i,j,I^{\prime },J^{\prime },s\right):\in \left\{a(i,j)\cap A(I^{\prime },J^{\prime })\right\}}{\sum e\left(i,j,I^{\prime },J^{\prime },s\right):\in \left\{A(I^{\prime },J^{\prime })\right\}}},\end{array}\end{array}$$

where the numerator is the sum of the emissions within the intersection of a(i,j) and $A(I^{\prime },J^{\prime })$ and the denominator is the sum of the emissions within $A(I^{\prime },J^{\prime })$. The resulting total concentration, using this method, may be higher or lower than the original EMEP concentrations because it reflects the impact of moving the EMEP grid in space. This is easiest to visualise if the moving window is the same size as the EMEP grid. If the moving window was centred on an area with concentrated emissions, that are in reality spread over two EMEP grids, then when using the emission weighting the new EMEP local contribution would be higher, the non-local lower and the total would be different; see Fig. 2. The opposite is also true if the moving window were placed over a region with low emissions, the local contribution would be lower and the non-local higher. Due to this, it is not possible simply to subtract the local EMEP contribution from the total to get the non-local EMEP contribution, as detailed in Eq. (5).

To address this, the non-local EMEP contribution is also calculated using the moving window with Eq. (9). The first term is the non-local contribution for a particular source and is calculated with the area weighting distribution since non-local contributions, those outside the local fraction region, do not have any associated emission or local fraction for weighting. An additional correction term, second term in Eq. (9), accounts for the non-local contributions from local contributions on the EMEP edge grids, those parts of the EMEP grids that are outside the moving window and not included as a local contribution in Eq. (6). In those cases, the local EMEP contribution outside the moving window must be converted to a non-local contribution and subtracted from the calculated non-local value (first term in Eq. 9).

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& C_{\mathrm{G},\mathrm{nonlocal}}\left(i,j,s\right)=\sum _{I^{\prime }=I-n_{\mathrm{mw}}/\mathrm{2}}^{I+n_{\mathrm{mw}}/\mathrm{2}}\sum _{J^{\prime }=J-n_{\mathrm{mw}}/\mathrm{2}}^{J+n_{\mathrm{mw}}/\mathrm{2}}\\ & C_{\mathrm{G},\mathrm{nonlocal}}\left(I^{\prime },J^{\prime },s\right)\cdot w_{\mathrm{a}}\left(i,j,I^{\prime },J^{\prime }\right)-\\ & \sum _{I^{\prime }=-n_{\mathrm{mw}}/\mathrm{2}\left(I^{\prime }\ne \mathrm{0}\right)}^{+n_{\mathrm{mw}}/\mathrm{2}}\sum _{J^{\prime }=-n_{\mathrm{mw}}/\mathrm{2}\left(J^{\prime }\ne \mathrm{0}\right)}^{+n_{\mathrm{mw}}/\mathrm{2}}\\ & \left(\begin{array}{l}{C}_{\mathrm{G},\mathrm{local}}\left(I,J,I^{\prime }\to I,J^{\prime }\to J,s\right)\cdot w\left(i,j,I^{\prime },J^{\prime },s\right)\\ +C_{\mathrm{G},\mathrm{local}}\left(I,J,I\to I^{\prime },J\to J^{\prime },s\right)\cdot w\left(i,j,I,J,s\right)\end{array}\right)\end{array}\end{array}$$

In Eq. (9), the weighting term w represents either the emission (w_e) or the area (w_a) weighting, depending on the choice of weighting method.

These local and non-local calculations are carried out for each emission source individually so the non-local contribution is also dependent on source and the non-local component for any particular source will also contain the local contributions from the other sources. This makes creating a final non-local contribution complicated. To solve this, all the source specific $C_{\mathrm{G},\mathrm{local}}+C_{\mathrm{G},\mathrm{nonlocal}}$ contributions are averaged and the sum of the C_G,local source contributions are subtracted to obtain the final C_G,nonlocal. The final non-local contribution at each subgrid C_SG,nonlocal (Eq. 1) is equivalent to the EMEP non-local C_G,nonlocal contribution and is calculated by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& C_{\mathrm{SG},\mathrm{nonlocal}}\left(i,j\right)={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{source}}}}\sum _{s=\mathrm{1}}^{n_{\mathrm{source}}}\left\{(C_{\mathrm{G},\mathrm{local}}\left(i,j,s\right)\right.\\ & \left.+C_{\mathrm{G},\mathrm{nonlocal}}\left(i,j,s\right))\right\}-\sum _{s=\mathrm{1}}^{n_{\mathrm{source}}}\left\{C_{\mathrm{G},\mathrm{local}}\left(i,j,s\right)\right\}.\end{array}\end{array}$$

In the case of area weighting, where the sum of local and non-local contribution is the same as the original EMEP total concentration, then the first term in the summation is equivalent to the original EMEP concentration without summation. The method is illustrated in one dimension in Fig. 2.

The calculation based on emission weighting is computationally more expensive than the area weighting and is only used when necessary and appropriate, e.g. when n_mw=1 and when subgrid and grid emissions are consistent with each other.

Figure 2 Illustration of the moving window interpolation method employed in uEMEP. Shown is the 1-D visualisation of the 2-D method described in Eqs. (6)–(10) for n_mw=1. For clarity, in the figure, the terms C_G,local and C_G,nonlocal are written as L and NL, respectively.

3 uEMEP model process description and parameterisation

In this section, uEMEP process parameterisations are described. In regard to the dispersion modelling, uEMEP is intended to integrate closely with EMEP. To enable this, dispersion schemes based on parameterisations used in EMEP have been implemented. In the Supplement, additional equations (Sects. S1–S3) are provided and a number of optional additional parameterisations are also described (Sect. S4).

3.1 Subgrid Gaussian dispersion modelling for hourly calculations

A standard Gaussian narrow plume dispersion model formulation, e.g. Seinfeld and Pandis (1998), is used in the subgrid dispersion calculations with multiple reflections from the surface (z=0) and boundary layer height (z=H). Generically, the Gaussian plume calculation can be written as

$$\begin{array}{}\text{(11)}& C\left(x,y,z\right)={\displaystyle \frac{Q}{U}}I(x,y,z),\end{array}$$

where for the sake of clarity we have dropped references to subgrid indexes as given in Sect. 2 and use coordinates instead of indices. Here, C is the concentration, Q is the emission source strength and I is the plume intensity given by

$$\begin{array}{}\text{(12)}& I(x,y,z)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_y{\mathit{\sigma }}_z}}\exp \left({\displaystyle \frac{-y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right)\sum _{i=\mathrm{1}}^{i=\mathrm{6}}\left\{\exp \left({\displaystyle \frac{-{\left(z-h_i\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)\right\}.\end{array}$$

Here, h_i represents the plume emission height and five additional virtual plume emission heights after single and double reflections from the surface and boundary layer top (H) given by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{h}_i=& \left[h_{\mathrm{emis}},-h_{\mathrm{emis}},\mathrm{2}H-h_{\mathrm{emis}},\mathrm{2}H+h_{\mathrm{emis}},-\mathrm{2}H\right.\\ & \left.+h_{\mathrm{emis}},-\mathrm{2}H-h_{\mathrm{emis}}\right].\end{array}\end{array}$$

For the well-mixed plume case, when σ_z is of the order of H, we define a threshold beyond which the plume concentration is constant throughout the boundary layer. This is specified to occur when σ_z>0.9H, leading to an intensity given by

$$\begin{array}{}\text{(14)}& I(x,y,z)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{y}H}}\exp \left({\displaystyle \frac{-y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right).\end{array}$$

The Gaussian dispersion parameters (σ_z and σ_y) used in Eq. (12) may be determined empirically (Smith, 1973; Martin, 1976; Turner, 1994; Liu et al., 2015) or through a range of methods based on theoretical and semi-empirical considerations (Seinfeld and Pandis, 1998). Venkatram (1996) also discusses the relationship between empirically and theoretically based dispersion parameters. Standard Gaussian plume models do not take into account variable vertical profiles of wind speed or diffusivity. Some non-Gaussian descriptions are available based on the application of power laws to these profiles and the vertical integration of the diffusion equation (Chaudhry and Meroney, 1973; van Ulden, 1978; Venkatram et al., 2013) but this then creates the problem of defining power laws that “fit” varying wind and dispersion profiles over the entire boundary layer. Instead of this, we use the centre of mass of the plume (z_cm) to define the height at which the advective wind speed and eddy diffusivity (K_z) are defined and allow this to vary dependent on the plume travel distance, giving a similar effect to the plume dispersion as the non-Gaussian vertically integrated derivation. A similar methodology is employed by the OPS model (Sauter et al., 2018). We then use a combination of eddy diffusivity (K_z) vertical profiles, Lagrangian timescales and centre of mass plume placement, along with initial values σ_z0 and σ_y0, to determine σ_z and σ_y values. The aim of this combination is to provide realistic plume dispersion over short distances but to asymptotically approach the same K_z values used in the EMEP model dispersion scheme over longer distances. In addition, the methodology is implementable at all emission heights and takes into account both surface roughness and atmospheric boundary layer height.

Following the methodologies outlined in Seinfeld and Pandis (1998), we describe the dispersion parameters σ_z and σ_y as a function of pollutant travel time from source (t) using

$$\begin{array}{}\text{(15a)}& {\displaystyle }{\mathit{\sigma }}_z\left(t\right)={\mathit{\sigma }}_{z\mathrm{0}}+\sqrt{\mathrm{2}{K}_z\left(z\right)t{f}_t}\text{(15b)}& {\displaystyle }{\mathit{\sigma }}_y\left(t\right)={\mathit{\sigma }}_{y\mathrm{0}}+\sqrt{\mathrm{2}{K}_y\left(z\right)t{f}_t},\end{array}$$

where σ_z0 and σ_y0 are initial dispersion parameters, K_z(z) and K_y(z) are the vertical profiles of vertical and horizontal eddy diffusivity, and f_t is a factor dependent on the Lagrangian integral timescale τ_l given by

$$\begin{array}{}\text{(16)}& f_t=\mathrm{1}+\left({\displaystyle \frac{{\mathit{\tau }}_l}{t}}\exp \left(-{\displaystyle \frac{t}{{\mathit{\tau }}_l}}\right)-\mathrm{1}\right).\end{array}$$

There are many varying methods for calculating the Lagrangian integral timescale (Seinfeld and Pandis, 1998; Hanna, 1981; Venkatram, 1984). We use the formulation from Hanna (1981):

$$\begin{array}{}\text{(17)}& {\mathit{\tau }}_l=\mathrm{0.6}{\displaystyle \frac{max\left(z_{\mathrm{emis}},z_{\mathit{\tau }\min }\right)}{u_{\ast }}},\end{array}$$

where z_emis is the emission height, u_∗ is the friction velocity, and $z_{\mathit{\tau }min}=\mathrm{2}$ m. Time is calculated from the advective velocity:

$$\begin{array}{}\text{(18)}& t={\displaystyle \frac{max\left(x_{min},x\right)}{U\left(z\right)}},\end{array}$$

where x_min is half a subgrid, U(z) is the vertical wind speed profile, and x is the down-plume distance.

In order to be compatible with the EMEP model, the same K_z vertical profile parameterisation is used in Eq. (15a) that is used in EMEP (Simpson et al., 2012). This parameterisation is provided in the Supplement (Eqs. S1–S2).

The centre of mass of the plume is calculated using the same Gaussian formulation with reflection as given in Eq. (12) by integrating the plume intensity over the boundary layer height (H) using

$$\begin{array}{}\text{(19)}& z_{\mathrm{cm}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{H}zI\left(z\right)\mathrm{d}z}{{\int }_{\mathrm{0}}^{H}I\left(z\right)\mathrm{d}z}}.\end{array}$$

This integral can be analytically solved to give

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{z}_{\mathrm{cm}}=& {\displaystyle \frac{{\mathit{\sigma }}_z}{\sqrt{\mathrm{2}\mathit{\pi }}}}\sum _{i=\mathrm{1}}^{i=\mathrm{6}}\exp \left({\displaystyle \frac{-h_i^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)-\exp \left({\displaystyle \frac{-{\left(H-h_i\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}}}\right)\\ & +{\displaystyle \frac{h_i}{\mathrm{2}}}\left(\mathrm{erf}\left({\displaystyle \frac{h_i}{\sqrt{\mathrm{2}}{\mathit{\sigma }}_z}}\right)+\mathrm{erf}\left({\displaystyle \frac{(H-h_i)}{\sqrt{\mathrm{2}}{\mathit{\sigma }}_z}}\right)\right),\end{array}\end{array}$$

and for the well-mixed case where

$$\begin{array}{}\text{(21)}& {\mathit{\sigma }}_z>\mathrm{0.9}H\text{ then }{z}_{\mathrm{cm}}=\mathrm{0.5}H.\end{array}$$

The vertical wind profile is calculated in a similar way to Gryning et al. (2007) based on decreasing turbulent shear with height.

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& U\left(z\right)={\displaystyle \frac{u_{\ast \mathrm{0}}}{\mathit{\kappa }}}\left(\log \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_m+\mathit{\kappa }{\displaystyle \frac{z}{z_l}}\left(\mathrm{1}-{\displaystyle \frac{z}{\mathrm{2}H}}\right)\right.\\ & \left.-{\displaystyle \frac{z}{H}}\left(\mathrm{1}+b{\displaystyle \frac{z}{\mathrm{2}L}}\right)\right)\text{for}L\ge \mathrm{0}\\ & U\left(z\right)={\displaystyle \frac{u_{\ast \mathrm{0}}}{\mathit{\kappa }}}\left(\log \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_m+\mathit{\kappa }{\displaystyle \frac{z}{z_l}}\left(\mathrm{1}-{\displaystyle \frac{z}{\mathrm{2}H}}\right)\right.\\ & \left.-{\displaystyle \frac{z}{H}}{\displaystyle \frac{\left(\left(az-L\right){\mathit{\varphi }}_m+L\right)}{a(p+\mathrm{1})}}\right)\text{for}L<\mathrm{0}\end{array}\end{array}$$

The stability functions (ψ_m and ϕ_m) are defined in Eqs. (S3)–(S4) in the Supplement, and the assumptions behind the wind profile derivation are given in Eqs. (S5)–(S8). There is no turning of the wind direction with height. Equation (22) is used to derive u_∗0, based on modelled 10 m wind speed, boundary layer height H, Monin–Obukhov length (L) and surface roughness length z₀. The vertical wind profile is then derived from this. A minimum wind speed of 0.5 m s⁻¹ for all dispersion calculations has been imposed.

The average of the plume centre of mass height at the receptor point and the emission height, $z_{\mathrm{av}}=\mathrm{0.5}(z_{\mathrm{cm}}+h_{\mathrm{emis}})$, is then used to determine the vertical diffusion K_z(z_av) as well as the wind speed U(z_av) for use in Eqs. (15) and (19). The entire set of equations (Eqs. 15–22) is solved iteratively to obtain the final σ_z value at the receptor point. This iteration converges swiftly and generally only two iterations are required.

The horizontal eddy diffusivity K_y is not determined in EMEP so an alternative is required. K_y can be classically related to K_z through the relationship

$$\begin{array}{}\text{(23)}& K_y\left(z\right)={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{v}}(z{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{w}}(z{)}^{\mathrm{2}}}}{K}_z\left(z\right),\end{array}$$

based on the concepts used to define K (Seinfeld and Pandis, 1998). Garratt (1994) provides expressions for the vertical profile σ_v and σ_w under unstable conditions where the ratio σ_v∕σ_w is around 1.85 in the surface layer but decreases to 1 in the mixed layer. Under stable conditions, Nieuwstadt (1984) provides local scaling where this ratio is close to 2. For the current application, we choose the ratio ${\mathit{\sigma }}_{\mathrm{v}}/{\mathit{\sigma }}_{\mathrm{w}}=\mathrm{2}$ and apply it over the whole boundary layer.

It is also possible within the modelling setup to use the simpler empirical formulations of σ_z and σ_y, as presented in Eq. (24) and shown in Table 1. This is useful for testing and comparison, and necessary when using the rotationally symmetric plume parameterisation, Sect. 3.2. See Seinfeld and Pandis (1998) for a presentation of these.

In Fig. 3, we show two example sets of σ_z curves for near-surface (1 m) and elevated (50 m) releases as calculated with the K_z methodology for three separate stabilities. For reference, the dispersion curves from ASME (American Society of Mechanical Engineers; Smith, 1973) are also shown. These often-used dispersion parameters are relevant for 1 h averaging times. The ASME σ_z curves are given in Pasquill stability classes and the conversion from their dependency on L and z₀ is achieved using the conversion methodology described by Golder (1972). Parameters used in the calculation of the three curves are provided in Table 1.

In addition to the parameterisations presented here, uEMEP also includes parameterisations, provided in the Supplement, for plume meandering and change of wind direction (Sect. S3.4.1), traffic-induced initial dispersion (Sect. S3.4.3) and road tunnel internal deposition and emissions (Sect. S3.4.5).

Figure 3 Comparison of derived σ_z curves discussed in the text with standard ASME curves (Smith, 1973) using Eq. (24). To the left is a 1 m release and to the right a 50 m release. Three different stability classes, specified by the Monin–Obukhov lengths (L), are shown. The K_z method is shown as a solid line and the ASME curves as dashed lines. The ASME curves have no release height or surface roughness dependence but are included as reference. Values of z₀=0.5 m, relevant for urban areas, and σ_z0=0 are used.

Table 1 Parameters used for calculating the curves shown in Fig. 3.

3.2 Rotationally symmetric Gaussian plume model for annual mean calculations

When applying uEMEP to annual mean emissions, a rotationally symmetric Gaussian plume is used. It is possible to derive an approximate analytical solution to the Gaussian plume equation assuming that wind directions are homogeneously distributed in all directions and that there is no strong dependence of wind speed or stability on wind direction. These conditions are usually not met but it is useful to have such a simplified analytical solution.

The starting point is the Gaussian plume model given in Eq. (12). In this case, we do not derive σ_y,z using the K_zvalue from EMEP but apply the commonly used power law formulation in order to derive an analytical solution:

$$\begin{array}{}\text{(24)}& {\mathit{\sigma }}_{(y,z)}={\mathit{\sigma }}_{\mathrm{0}(y,z)}+a_{(y,z)}{x}^{b_{(y,z)}}.\end{array}$$

Values for the dispersion parameters a and b may be taken from the literature (Seinfeld and Pandis, 1998) but we use the ASME curves (Smith, 1973) under neutral conditions to specify these.

The rotationally symmetric version of this equation can be derived by rewriting the equation in cylindrical coordinates with appropriate approximations (second order), based on the slender plume assumption, and integrating over all angles. The resulting rotationally symmetric intensity I_rot(r,z) as a function of r and z is then written

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& I_{\mathrm{rot}}\left(rz\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }\sqrt{\mathrm{2}\mathit{\pi }}r{\mathit{\epsilon }}_z\sqrt{\mathrm{1}+B}}}\mathrm{erf}\left({\displaystyle \frac{\mathit{\pi }\sqrt{\mathrm{1}+B}}{\sqrt{\mathrm{2}}{\mathit{\epsilon }}_{\mathit{\theta }}}}\right)\\ & \sum _{i=\mathrm{1}}^{i=\mathrm{6}}\left\{\exp \left({\displaystyle \frac{-{\left(z-h_i\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\epsilon }}_z^{\mathrm{2}}}}\right)\right\},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(26a)}& {\displaystyle }{\mathit{\epsilon }}_z={\mathit{\sigma }}_{\mathrm{0}z}+a_z{r}^{b_z}\text{(26b)}& {\displaystyle }{\mathit{\epsilon }}_{\mathit{\theta }}={\displaystyle \frac{\mathrm{1}}{r}}\left({\mathit{\sigma }}_{\mathrm{0}y}+a_y{r}^{b_y}\right)\text{(26c)}& {\displaystyle }B=-{\mathit{\epsilon }}_{\mathit{\theta }}^{\mathrm{2}}\left({\displaystyle \frac{b_z\left({\mathit{\epsilon }}_z-{\mathit{\sigma }}_{\mathrm{0}z}\right)}{r{\mathit{\epsilon }}_{\mathit{\theta }}}}+{\displaystyle \frac{b_y\left(r{\mathit{\epsilon }}_{\mathit{\theta }}-{\mathit{\sigma }}_{\mathrm{0}y}\right)}{{\mathit{\epsilon }}_z}}\right).\end{array}$$

The term B can be less than −1, typically when $r<\mathrm{2}{\mathit{\sigma }}_{\mathrm{0},y}$, which can lead to imaginary solutions. This is due to the second-order approximation made in converting to cylindrical coordinates. In that case, we write a second-order approximation based on Taylor series expansion around $B=-\mathrm{1}$ as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& I_{\mathrm{rot}}\left(r,z\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }r{\mathit{\epsilon }}_z{\mathit{\epsilon }}_{\mathit{\theta }}}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\pi }}^{\mathrm{2}}(\mathrm{1}+B)}{\mathrm{6}{\mathit{\epsilon }}_{\mathit{\theta }}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\pi }}^{\mathrm{4}}(\mathrm{1}+B{)}^{\mathrm{2}}}{\mathrm{40}{\mathit{\epsilon }}_{\mathit{\theta }}^{\mathrm{4}}}}\right)\\ & \sum _{i=\mathrm{1}}^{i=\mathrm{6}}\left\{\exp \left({\displaystyle \frac{-{\left(z-h_i\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\epsilon }}_z^{\mathrm{2}}}}\right)\right\}\text{ for}B<-\mathrm{1}.\end{array}\end{array}$$

A similar derivation has been carried out by Green (1980) using different assumptions for the form of Eq. (24).

3.3 Initial dispersion

In Sect. 3.1 and 3.2, the hourly and annual dispersion parameterisations are described. In both cases, initial values for σ_0(y,z) are required. Since we treat the sources as small area emitters, we set the initial σ_0y to correspond to these areas. A value of ${\mathit{\sigma }}_{\mathrm{0}y}=\mathrm{\Delta }y/\sqrt{\mathrm{2}\mathit{\pi }}\approx \mathrm{0.8}\left(\mathrm{\Delta }y/\mathrm{2}\right)$ will give a maximum subgrid centre concentration equivalent to the concentration that would be found if the emissions were distributed evenly in the subgrid. We then write the total initial dispersion to be

$$\begin{array}{}\text{(28)}& {\mathit{\sigma }}_{\mathrm{0}y}={\mathit{\sigma }}_{\mathrm{init},y}+\mathrm{0.8}{\displaystyle \frac{\mathrm{\Delta }y}{\mathrm{2}}}.\end{array}$$

In all applications of uEMEP, Δx=Δy. The other parameter, σ_init,y, is a specific initial dispersion width for each individual emission source, for example, 2 m for traffic and 5 m for shipping, heating and industry. This is generally much smaller than the emissions grids.

The initial value for σ_0z is also a combination of a specific emission initial dispersion, for example, ${\mathit{\sigma }}_{\mathrm{init},z}=\mathrm{5}$ m for residential wood combustion but also uses the displacement technique for the plume where the start of the plume is displaced upwind by Δx∕2 allowing the plume to grow vertically over half the subgrid distance. Tunnel exits are given an initial ${\mathit{\sigma }}_{\mathrm{init},z}=\mathrm{6}$ m to represent the extended size of the tunnel portals.

3.4 NO₂ chemistry for hourly means

The only chemistry included in uEMEP is the NO_x, O₃ chemical reactions. We use a similar methodology to Benson (1984, 1992) known as the discrete parcel method but use a weighted timescale over which the reactions take place. The following chemical reactions are involved, with O_x (O₃+NO₂) and NO_x (NO+NO₂) concentrations being conserved:

$$\begin{array}{}\text{(29a)}& {\displaystyle }\mathrm{NO}+{\mathrm{O}}_{\mathrm{3}}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}\text{(29b)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}+h\mathit{\nu }\to \mathrm{NO}+\mathrm{O}\text{(29c)}& {\displaystyle }{\mathrm{O}}_{\mathrm{2}}+\mathrm{O}+M\to {\mathrm{O}}_{\mathrm{3}}+M.\end{array}$$

Equation (29c) occurs on timescales much faster than the two other reactions and is taken to be instantaneous. The differential equation for the concentration of [NO₂] as a function of time is written as

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NO}}_{\mathrm{2}}\right]}{\mathrm{d}t}}=k_{\mathrm{1}}\left[\mathrm{NO}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right]-J\left[{\mathrm{NO}}_{\mathrm{2}}\right],\end{array}$$

where the concentrations are expressed in terms of molecules per cm³ and J is the photolysis rate (s⁻¹) for Eq. (29b) taken from the EMEP model (Simpson et al., 2012). The reaction rate k₁ for Eq. (29a) is given by

$$\begin{array}{}\text{(31)}& k_{\mathrm{1}}=\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{12}}\exp \left({\displaystyle \frac{-\mathrm{1310}}{T_{\mathrm{air}}}}\right)\left({\mathrm{cm}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right),\end{array}$$

as in the EMEP model and where T_air is in the atmospheric temperature (K).

We rewrite Eq. (30) in terms of the dimensionless ratios:

$$\begin{array}{}\text{(32a)}& {\displaystyle }{f}_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}{\left[{\mathrm{NO}}_x\right]}}\text{and}{f}_{{\mathrm{O}}_x}={\displaystyle \frac{\left[{\mathrm{O}}_x\right]}{\left[{\mathrm{NO}}_x\right]}}\text{(32b)}& {\displaystyle }{J}^{\prime }={\displaystyle \frac{J}{k_{\mathrm{1}}\left[{\mathrm{NO}}_x\right]}}\text{(32c)}& {\displaystyle }{t}^{\prime }=t{k}_{\mathrm{1}}\left[{\mathrm{NO}}_x\right].\end{array}$$

Equation (30) then becomes

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\mathrm{d}{f}_{{\mathrm{NO}}_{\mathrm{2}}}}{\mathrm{d}{t}^{\prime }}}=\left(\mathrm{1}-f_{{\mathrm{NO}}_{\mathrm{2}}}\right)\left(f_{{\mathrm{O}}_x}-f_{{\mathrm{NO}}_{\mathrm{2}}}\right)-J^{\prime }{f}_{{\mathrm{NO}}_{\mathrm{2}}}.\end{array}$$

The solution to this equation is

$$\begin{array}{}\text{(34)}& f_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{B}{\mathrm{2}}}{\displaystyle \frac{\left(\mathrm{1}-A\exp \left(B{t}^{\prime }\right)\right)}{\left(\mathrm{1}+A\exp \left(B{t}^{\prime }\right)\right)}}+{\displaystyle \frac{C}{\mathrm{2}}},\end{array}$$

where

$$\begin{array}{}\text{(35a)}& {\displaystyle }A={\displaystyle \frac{B+C-\mathrm{2}{f}_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{0}}}{B-C+\mathrm{2}{f}_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{0}}}}\text{(35b)}& {\displaystyle }B=\sqrt{C^{\mathrm{2}}-\mathrm{4}{f}_{{\mathrm{O}}_x}}\text{(35c)}& {\displaystyle }C=\mathrm{1}+f_{{\mathrm{O}}_x}+J^{\prime },\end{array}$$

and $f_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{0}}$ is the initial NO₂ fraction at time $t^{\prime }=\mathrm{0}$.

This solution is valid for a box model without dilution through dispersion since it does not take into account how changing NO_x and O_x concentrations over the plume travel time will affect the reaction rates. Though this could be accounted for when applied to a single source with assumed dilution rates, by adding a time-dependent diluting term to Eq. (30), this is not practically possible for multiple sources of differing dilutions. The concentrations of NO₂ at the start of the plume will be correctly calculated but NO₂ concentrations further from the plume will be slightly underestimated, since they do not have the higher initial reaction rates. Eventually, the concentrations will reach photostationary equilibrium and here too NO₂ will be correctly calculated. This special case for photostationary equilibrium in Eq. (35) occurs when $t^{\prime }\to \mathrm{\infty }$ and Eq. (34) becomes

$$\begin{array}{}\text{(36)}& f_{{\mathrm{NO}}_{\mathrm{2}}}={\displaystyle \frac{C-B}{\mathrm{2}}}.\end{array}$$

The non-linear nature of Eq. (34) also means that it cannot be consistently applied to Gaussian models since the shape of the plume will change due to the non-linearity. Despite this, this formulation is more physically realistic than the photostationary assumption often used in local-scale air quality modelling or other less physical parameterisations based on empirical fits. See Denby (2011) for an overview of the various NO₂ chemistry parameterisation methods used with Gaussian modelling. In order to calculate Eq. (34) in the model application, an initial NO_x and O_x concentration must be used and a travel time defined. For multiple sources, this travel time will vary so for each calculated subgrid concentration of NO_x from each contributing subgrid source (n_s sources) a travel time, t_s, is calculated based on the distance and wind speed. This is weighted based on the contribution of each source to the total subgrid NO_x concentration. This provides a final weighted travel time t_w that is applied in Eq. (34). This ensures that the nearest of the contributing subgrids, often with the highest contributing NO_x concentrations, are given a higher travel time weight. A minimum distance, and hence time, of half a subgrid is applied when calculating travel times.

$$\begin{array}{}\text{(37)}& t_{\mathrm{w}}={\displaystyle \frac{\sum _{s=\mathrm{1}}^{n_{\mathrm{s}}}{t}_{\mathrm{s}}{\left[{\mathrm{NO}}_x\right]}_{\mathrm{s}}}{\sum _{s=\mathrm{1}}^{n_{\mathrm{s}}}{\left[{\mathrm{NO}}_x\right]}_{\mathrm{s}}}}\end{array}$$

Comparisons with EMEP NO₂ calculations show that this chemistry scheme matches the results obtained by EMEP over longer time periods.

3.5 NO₂–NO_x conversion for annual means

When annual mean data are used, the hourly mean formulation cannot be applied. Instead, we use an empirically based conversion of NO_x to NO₂ based on the type of formulation from Romberg (1996) and updated by Bächlin and Bösinger (2008). A total of 3 years of Norwegian NO₂ measurements, 82 measurements in all, have been used to determine this relationship (Fig. 4).

$$\begin{array}{}\text{(38)}& \left[{\mathrm{NO}}_{\mathrm{2}}\right]={\displaystyle \frac{a\left[{\mathrm{NO}}_x\right]}{\left[{\mathrm{NO}}_x\right]+b}}+c\left[{\mathrm{NO}}_x\right]\end{array}$$

The fitted constants are determined to be a=20, b=30, and c=0.23. The estimated uncertainty in this conversion is around 10 %, based on the normalised root mean square error of the fitted and observed NO₂ concentrations.

Figure 4 NO₂ versus NO_x annual mean concentrations for all stations in Norway in the period 2013–2015. The fitted curve is based on Eq. (38).

This empirical relationship will vary from region to region, largely due to differences in O₃ concentrations and photolysis rates that are not included as part of the parameterisation. If used over large regions, for example, Europe, then the uncertainty in the NO₂ conversion will increase.

3.6 Implementation

3.6.1 Subgrid domains

Within uEMEP, individual domains are defined with differing resolutions and sizes, dependent on which modelling parameter is represented. Separate domains and subgrid resolutions are defined for each of the emission sources, for the time profiles of each emission source, for the meteorological data, for the population data and for the receptor subgrid concentrations. None of these are required to have the same resolution or size; however, the highest-resolution emission subgrid will define the receptor subgrid resolution, since there is no need to calculate on higher-resolution subgrids than is provided by the emissions. For emission subgrids with lower resolution than the final receptor subgrid domain, the dispersion calculations are first carried out at the same resolution as the emission domain and then bilinearly interpolated to the receptor subgrid. For most urban applications, this means that the choice of traffic subgrid resolution defines the highest-resolution subgrid.

Emission subgrids also contain properties for the dispersion calculations, such as initial dispersion parameters and emission heights. Each emission subgrid has only one emission height h_emis, one σ_init,z0 and one σ_init,y0 for each emission source type. When multiple emissions from the same source type are placed in an emission subgrid, the emission parameters are weighted by each individual emission. This is most relevant for industrial emissions which may have different emission heights from separate stacks within a single emission subgrid.

3.6.2 Selective subgrid calculations

uEMEP does not necessarily calculate concentrations at all receptor subgrids. Only subgrids which are within 3σ_y of a plume centre line will be calculated and also downwind selection is used (Sect. S3.4.2 in the Supplement). In addition, a number of selections can be made, allowing quicker calculations for particular applications. These include the following:

1. For calculation at defined receptor points, usually corresponding to measurement stations, uEMEP calculates the surrounding nine subgrids and uses bilinear interpolation to extract the concentrations at the required receptor position.

2. For calculation at population grids, concentrations will only be calculated at grids with non-zero population. This provides quicker exposure calculations than if the entire region was calculated.

3. A routine for selecting a higher density of subgrids near sources may also be used to speed up calculations. This applies most often to traffic emission subgrids that are near the surface and with large gradients near the source. This is less useful for higher release sources as their maximum impact occurs further downwind than their emissions. After calculation, the lower-density receptor subgrids are interpolated into the rest of the receptor subgrids, providing a full receptor subgrid domain.

3.6.3 Model inputs and outputs

Input data come from a variety of sources and the formatting of these sources varies. Emission input data are generally in text format whilst meteorological files are read from NetCDF files.

Output of the model is in the form of NetCDF files for either gridded data or point data, if receptor points have been defined. In both of these files, output includes the total concentrations of the pollutants along with the source contribution from each of the emission sources used in the calculation. Speciation of PM from EMEP can also be included in the output files, along with emissions, meteorology and population data.

4 Implementation in the Norwegian air quality forecast and analysis system

Though uEMEP has been applied in a number of applications we select the Norwegian forecast and analysis system (Norwegian air quality forecasting service, 2020) as an example. This application started operationally in the winter of 2018–2019 and provides daily forecasts of air quality for all of Norway 2 d in advance at subgrid resolutions of between 250 and 50 m. In addition, the same system is used to calculate air quality retrospectively for analysis and planning applications (Norwegian air quality expert user service, 2020). The compounds PM_2.5, PM₁₀, NO₂, NO_x and O₃ are calculated. For each of these, the local source contributions are determined separately for traffic exhaust, traffic non-exhaust, residential wood combustion, shipping and industry. A cascade of models is used, starting with EMEP MSC-W at 0.1^∘ European domain, EMEP MSC-W at 2.5 km Scandinavian domain and uEMEP 250–50 m Norway only (Fig. 5).

Figure 5 Model domain for the European and Scandinavian EMEP MSC-W calculations and the uEMEP calculations (© kartverket/norgeskart.no).

4.1 Calculation steps

We describe below the implementation steps used in the Norwegian forecasting and analysis system. This implementation of uEMEP uses the independent emission and replacement downscaling method (method 2 in Sect. 2). The following steps are carried out:

• 1.

  High-resolution emission data for Norway are calculated for each forecast (Sect. 4.2) and are aggregated into the EMEP MSC-W Scandinavian model grid. Some of these emissions require meteorological data.

• 2.

  The EMEP MSC-W model is used to calculate the large-scale concentration distribution on an hourly basis, nesting from a European domain (∼ 0.1^∘) to the Scandinavian domain (2.5 km) (Fig. 5). Within Norway, the aggregated high-resolution aggregated emissions are implemented. Both EMEP calculations provide the local fraction (Sect. 2.3) in a region of 5×5 EMEP grids.

The following three steps are then undertaken to calculate the uEMEP concentrations:

• 3.

  For the Norwegian forecast system, the entire country is split into 1864 separate tiles of varying sizes and resolutions, with the resolution depending on the population and emission sources within each tile. Tiles with resolutions of 250 m can be as large as 40×40 km², whilst tiles with resolutions of 50 m are no larger than 5×5 km². Tiling the calculations is a form of external parallelisation and is optimised for both runtime and memory use. A 2 d forecast run on 196 processors takes roughly 1 h of CPU time.

• 4.

  The high-resolution emission data from the various source sectors (Sect. 4.2) are placed into the emission subgrids in uEMEP. These are between 50–250 m in width, depending on the emissions available and on the population density of the region. Emission grid domains extend beyond the size of each tile so that the calculations are consistent over tile borders.

• 5.

  uEMEP Gaussian dispersion modelling is applied (Sect. 3.1) using the subgrid emissions as sources and the concentrations are calculated at each subgrid. Only subgrid emissions within a region defined by a 4×4 EMEP grid area are included in the subgrid calculation, i.e. 10×10 km², corresponding to the extent of the moving window. This 4×4 limit guarantees that the calculation will always be carried out within the EMEP 5×5 local fraction region.

The final steps combine the EMEP gridded concentrations with the uEMEP subgrid concentrations in the following way:

• 6.

  At each subgrid, the non-local contribution from the neighbouring 4×4 EMEP grids is calculated (Sect. 2.4). The calculation is carried out for each source sector and each primary compound.

• 7.

  The uEMEP calculations are then added to the non-local EMEP concentrations. In the case of PM, all non-primary species are also added to the local and non-local EMEP primary concentrations.

• 8.

  For NO₂, the chemistry (Sect. 3.4) is applied to determine NO₂ and ozone for each subgrid.

• 9.

  Subgrid concentrations and their contributions are saved along with the PM speciation from EMEP in NetCDF format.

• 10.

  The forecasts are made available to a public website through an API and Web Map Tile Server (Norwegian air quality forecasting service, 2020).

The system is schematically illustrated in Fig. 6. The following sections describe some steps in more detail.

Figure 6 Flow diagram showing the various components of the Norwegian EMEP/uEMEP forecast system.

4.2 Emissions

The EMEP calculations make use of the CAMS-REG-AP_v1.1 regional anthropogenic emission dataset everywhere in Europe (Kuenen et al., 2014; Granier et al., 2019). Only in the 2.5 km Scandinavian calculation, and only in Norway, are the emissions replaced with the aggregated high-resolution dataset. The alternative emissions used in the calculations for Norway are

• road traffic exhaust emissions,

• road traffic non-exhaust emissions,

• residential wood combustion,

• shipping and

• industry.

These emission sources are described in Sect. S4.2 in the Supplement. For other sectors, the CAMS-REG-AP_v1.1 emissions are also used in Norway, but these emissions are not downscaled using uEMEP.

4.3 Meteorology

The meteorological forecast data used for the European EMEP model calculations are based on the Integrated Forecasting System (IFS, 2020) from the European Centre for Medium-Range Weather Forecasts (ECMWF, 2020). The Scandinavian EMEP model calculation uses the AROME-MetCoOp model for modelling meteorology over Scandinavia (Müller et al., 2017). This last model calculates meteorology at a resolution of 2.5 km and provides forecasts for 66 h in advance. The EMEP MSC-W Scandinavian domain uses the same gridding and projection as the meteorological forecast model but in a smaller domain.

4.4 EMEP MSC-W model implementation

The European EMEP MSC-W model calculation is based on the same daily forecast provided for the Copernicus Atmosphere Monitoring Service (CAMS, 2020; Tarrason et al., 2018) but is run independently and provides boundary conditions for the Scandinavian implementation of EMEP MSC-W at 2.5 km. The Scandinavian EMEP MSC-W calculation includes the Norwegian emission sources described in Sect. 4.2 and also delivers the necessary local fraction information for use in uEMEP.

4.5 uEMEP model implementation

uEMEP calculates concentrations for all of Norway on grids with resolutions between 50 and 250 m using 1864 individual tiles as described in Sect. 4.1. The resolution of these tiles is defined by the population density and road density information. Tiles with higher population density use 50 m resolution, whilst tiles with lower population density but some traffic have a resolution of 125 m. Tiles with very low traffic density but with shipping or wood burning emissions have a resolution of 250 m, corresponding to the emission resolution. Separate calculations are carried out at measurement sites, 72, with a subgrid resolution of 25 m. An example of a PM₁₀ forecast is shown in Fig. 7.

Figure 7 Example maps of PM₁₀ concentrations taken from the forecast on 24 February 2020 at 18:00 UTC. The resolution in populated city regions is 50 m. High PM₁₀ concentrations along roads are mainly the result of road dust emissions (© kartverket/norgeskart.no).

5 Results

5.1 Validation against observations for the Norwegian forecasting and assessment system

Here, we present a visual summary of results for NO₂, PM₁₀ and PM_2.5 for the year 2017, when there were 72 operational air quality stations. Not all stations measure all components so the total number of available stations for NO₂ and PM with coverage of more than 75 % is between 34 and 45. The station positions are shown in Fig. 8.

Figure 8 Positions of all 72 monitoring stations in Norway: 33 for PM_2.5, 36 for NO₂, 45 for PM₁₀ and 8 for O₃ (© kartverket/norgeskart.no).

Figure 9 shows the comparison of modelled and observed NO₂ for annual average at each station (scatter plot) and daily mean temporal profile averaged over all stations. Included in the scatter plot are the Scandinavian EMEP MSC-W results at 2.5 km. The spatial correlation is quite high, r²=0.81 for uEMEP with little negative bias (FB of −5.9 %). The temporal variation over the whole year is also well represented when averaged over all stations (r²=0.79).

Figure 9 (a) Scatter plot comparison of modelled and observed NO₂ for annual average concentrations at each station. Shown are the results for both uEMEP and EMEP. Comparative statistics are shown for the uEMEP calculation. EMEP and observed means are also included. (b) Daily mean temporal profile averaged over all stations. Source contributions are shown for the temporal modelling results along with the EMEP 2.5 km calculation (EMEP4NO). In total, 36 stations are used in the comparison.

Figure 10 shows the comparison of modelled and observed PM₁₀ for annual average at each station (scatter plot) and daily mean temporal profile averaged over all stations. Included in the scatter plot are the Scandinavian EMEP results at 2.5 km. The spatial correlation is low, r²=0.29 for uEMEP with little negative bias (FB of −9.2 %). The temporal variation over the whole year is well represented when averaged over all stations (r²=0.61) but the model has a negative bias of 4 µg m³ over much of the summer period. Road dust events in the springtime are well captured by the NORTRIP emission model.

Figure 10 (a) Scatter plot comparison of modelled and observed PM₁₀ for annual average concentrations at each station. Shown are the results for both uEMEP and EMEP. Comparative statistics are shown for the uEMEP calculation. EMEP and observed means are also included. (b) Source contributions from both uEMEP and EMEP models are shown for the temporal modelling results along with the EMEP 2.5 km calculation (EMEP4NO). A total of 45 stations are used in the comparison.

Figure 11 shows the comparison of modelled and observed PM_2.5 for annual average at each station (scatter plot) and daily mean temporal profile averaged over all stations. Included in the scatter plot are the Scandinavian EMEP results at 2.5 km. The spatial correlation is good, r²=0.49 for uEMEP with little negative bias (FB $=-\mathrm{10.5}$ %). The temporal variation over the whole year is well represented when averaged over all stations (r²=0.67) but the model has a negative bias of 2 µg m³ over much of the summer period. Residential wood combustion (heating) events in the winter are well captured by the emission model MetVed.

Figure 11 (a) Scatter plot comparison of modelled and observed PM_2.5 for annual average concentrations at each station. Shown are the results for both uEMEP and EMEP. Comparative statistics are shown for the uEMEP calculation. EMEP and observed means are also included. (b) Source contributions are shown from both uEMEP and EMEP models for the temporal modelling results along with the EMEP 2.5 km calculation (EMEP4NO). A total of 33 stations are used in the comparison.

5.2 Verification and sensitivity tests

In addition to the validation against monitoring data, a number of verification and sensitivity experiments have been undertaken with the model. These include

• comparison of single annual mean calculations with the mean of hourly calculations,

• sensitivity to the moving window size,

• sensitivity to the choice of resolution,

• sensitivity to the temperature dependence of NO_x exhaust emissions and

• sensitivity to the choice NO₂ ∕ NO_x initial exhaust ratio.

These sensitivity tests are provided in the Supplement (Sect. S5). We present only conclusions from these.

5.2.1 Comparison of single annual mean calculations with the mean of hourly calculations

In Sect. 3, we describe two methods for calculating dispersion. One is based on the hourly meteorological and emission data (Sect. 3.1) and the other on annual mean data (Sect. 3.2). A rotationally symmetric dispersion kernel (Eq. 25) is used for dispersion of the annual mean emissions. Tests using the same dispersion parameters in both annual and hourly calculations (Sect. S5.1) give very similar results for both methodologies indicating the validity of the annual mean calculation. When K_z-based dispersion (Eqs. 15–23) is used in the hourly calculations, there is a larger difference between the two methods because of the difference between the two dispersion parameterisations. We conclude that the time-saving advantage of the single annual mean calculation, approximately 10 000 times faster, and the similarity to the hourly calculation make it an efficient and valid method for calculating high-resolution annual maps of air quality.

5.2.2 Sensitivity to the moving window size

The size of the moving window region within which uEMEP calculates local high-resolution concentrations should impact on the results since smaller moving windows will include less locally modelled contributions and more non-local EMEP contributions. This has been verified in a sensitivity study, (see Sect. S5.2 in the Supplement). In this sensitivity experiment, the moving window size was varied from n_mw=1 to eight EMEP grids, and calculations were made at existing measurement sites. The mean concentrations are shown to be quite insensitive to the choice of this region, particularly for PM₁₀. Generally, the reduction in the local contribution is well balanced with the increase in the non-local contribution when reducing the size of the moving window, verifying the methodology. It is recommended to use a minimum of two EMEP grids for the moving window region.

5.2.3 Sensitivity to the choice of resolution

The choice of subgrid resolution will impact on the calculated concentrations, both in concentration levels and in spatial distribution. An experiment where a range of subgrid resolutions were tested, from 15 to 250 m, was carried out (Sect. S5.3). Calculations were made at the positions of the Norwegian monitoring sites, most of which are traffic sites. The results showed that even at resolutions of 250 m the mean concentrations for all stations were very similar. At 100 m resolution, compared to the reference of 25 m, the difference in annual mean was no larger than 15 % at any one station with a normalised root mean square error (NRMSE) of 6 %. The NRMSE increased to 11 % for the 250 m calculation with a maximum deviation of 40 % at one station. We conclude that 100 m resolution will provide good concentration estimates for near-road calculations though higher resolutions may be required, depending on the application.

5.2.4 Sensitivity to the temperature dependence of NO_x exhaust emissions

The temperature dependence of the NO_x traffic exhaust emissions was assessed by running the model with and without this dependency (Sect. S5.4). With this correction, the results show a significant improvement in the station mean time series correlation (from r²=0.68 to 0.79) and improved correlation in both the daily (from r²=0.56 to 0.60) and annual (from r²=0.76 to 0.78) mean calculations. Bias is also reduced from −20 % to −3 %. The correction factor used (Eq. S13) still requires further evaluation and should be considered only as an initial estimate.

5.2.5 Sensitivity to the choice NO₂ ∕ NO_x initial exhaust ratio

In the calculations shown in Fig. 9 for NO₂, an initial NO₂ ∕ NO_x exhaust emission ratio of 0.25 was used. This reflects the large portion of diesel vehicles used in Norway and the high NO₂ ∕ NO_x ratio of these (Hagman et al., 2011). However, comparison of modelled ratios of NO₂ ∕ NO_x indicate this ratio may be too high. This was assessed by running the model with three different ratios: 0.15, 0.25 and 0.35. The results (Sect. S5.5) show that an NO₂ ∕ NO_x ratio of 0.15 most closely fits the observed ratio and this ratio will be implemented in further applications of the model.

6 Discussion

The aim of uEMEP is to provide downscaling capabilities for the EMEP MSC-W model with the intention of improving exposure estimates and more realistic concentrations at high resolution over large areas. The example application provided, the Norwegian air quality forecast and expert user service, is an example of how high-resolution coverage over large regions (countries) can be achieved. The validation carried out in Sect. 5.1 shows that the modelling system provides moderate to good comparison with observations. The best results are for NO₂, chiefly because we have the best information concerning emissions that contribute to these concentrations, i.e. traffic exhaust. The lower correlation of PM is indicative of the difficulties in modelling emissions such as residential wood burning and road dust emissions. That NO₂ is well modelled indicates that the problems lie largely with the emissions rather than the dispersion model itself. In addition, a large proportion of PM is due to medium to long-range transport and secondary formation of particles. This is not part of the uEMEP model but relies on the EMEP MSC-W model and the emissions included there.

The strength of the modelling system is in the integration of uEMEP with EMEP through the use of the local fraction. This allows downscaling anywhere within an EMEP domain provided that suitable proxy data are available for the downscaling. This is an important aspect of the modelling and is the link that can bind the regional- and local-scale emission communities. Usually the proxies used for regional-scale emission inventories are not available to the user so that exactly how these emissions are made, quantitatively, is unknown to the user. In addition, as the resolution of regional-scale emission inventories increases, so too does the need for improved spatial distribution proxies. Population density, successfully used to distribute a range of emission sectors on low-resolution grids (> 10 km), is no longer suitable for many sectors since at high resolution the emissions are no longer correlated with population. This was discussed in an earlier paper (Denby et al. 2011) and remains problematic.

When implementing uEMEP, it is highly desirable that the emissions used in both uEMEP and EMEP models are consistent with one another. This has been achieved for the Norwegian application for the traffic, domestic heating, shipping and industry sectors. However, other sources, such as other mobile combustion sources associated with construction and other activities, are not included. These can be of importance locally even if they are not significant on the regional scale. There is no clear methodology available on how to implement these emissions at the required resolution.

The modelling system has limitations. Currently, only primary emissions, with the exception of NO₂ formation, are dealt with. Some secondary formation of particles will likely occur within the local region used for the uEMEP model domain but these are not currently accounted for. uEMEP is also a Gaussian model that does not take into account obstacles of any type. When achieving resolutions of 50 m, buildings start to play an important role in the transport and dispersion. The region covered by the local-scale modelling, the moving window region, is necessarily limited in extent. Sensitivity studies show that this has limited impact on mean concentrations but for source sectors such as industry, that are released at height, the limited calculation region may not be large enough to include all of the plumes impact region.

In many ways, the increase in resolution to almost street level puts new demands on the modelling system that were not necessary to consider previously. For regional-scale modellers, the downscaling can provide considerable improvement to regional calculations. However, from a local-scale modelling perspective, the local-scale information may not be of sufficient quality to be useful to local users. This is most important when only proxy data are available for downscaling rather than actual bottom up emissions. In the end, if high-resolution modelling is to be used at the local scale, then similar high-quality emission data will be required if the results are to be useful to users.

There are a number of aspects of the modelling system that can be, and are being, improved. These include

• implementation of dry and wet deposition in uEMEP (currently not included in this version),

• improving the annual mean dispersion kernel dispersion parameters to be more consistent with the hourly K_z methodology,

• implementing necessary secondary formation of PM in uEMEP,

• further assessment of the K_z Gaussian plume methodology and

• refinement of the temperature dependence of NO_x traffic exhaust emissions.

A number of aspects were not treated in this paper but will be topics of further studies. These include population exposure and the impact of resolution, trend assessment in emissions and analysis of road dust emissions for all of Norway. In addition, the modelling system is being applied in a number of different countries and results of these applications will be further described and assessed.

7 Conclusions

This paper presents and documents a new downscaling model and method (uEMEP) for use in combination with the EMEP MSC-W chemical transport model. Process descriptions and parameterisations within the uEMEP model are provided and the methodology for combining uEMEP with EMEP MSC-W local fraction calculations is elaborated. An example application, the Norwegian air quality forecast system, is presented and validation for the year 2017 at all available Norwegian air quality stations is provided. A number of verification and sensitivity studies are summarised in the paper and expanded in the Supplement.

The uEMEP model provides a new methodology for downscaling regional-scale chemical transport models but can currently only be applied together with the EMEP MSC-W model since this is the only model with the necessary local fraction calculation. The uEMEP model is based on Gaussian modelling that has existed for many years but it does use specific parameterisations to describe the dispersion parameters in order to be compatible with the EMEP model application.

uEMEP can provide improved exposure estimates if suitable proxy data for emissions are available and can be applied to regions as large as the regional-scale CTM in which it is imbedded. It can also represent concentrations down to street level, though not street canyons, and is comparable with traffic monitoring sites. This makes it a unique system for assessment, policy application and forecasting purposes.