Model Description

This page includes a description of the uEMEP model.

Introduction

Methodology

uEMEP can use two downscaling methods, both of which are employing Gaussian dispersion modeling to achieve redistribution of concentrations in high resolution (Denby et al.[1]). Both methods also utilize the “local fraction” feature in the EMEP model (Wind et al.[2]) to avoid double counting of emissions and to allow near-seamless integration of the two models:

1. Emission redistribution

2. Independent emissions

Emission redistribution is used when only approximate proxy data for emissions are available. In this case, emissions from EMEP are redistributed based on proxies such as population density, road network or landuse data. Independent emissions are used when high resolution emission data are available which are consistent between EMEP and uEMEP.

Subgrid calculation method

The total subgrid contrations, denoted as \(C_{SG}(i,j)\) in subgrid \((i,j)\) are calculated by summing the local component, \(C_{SG,local}(i,j)\), and the non-local component, \(C_{SG,nonlocal}\), derived from the EMEP grid concentration:

(1)\[C_{SG}(i,j) = C_{SG,local}(i,j) + C_{SG,nonlocal}(i,j)\]

Here, \(SG\) represents the uEMEP subgrid value, while \(G\) denotes EMEP grid values in the subsequent equations.

The local part of the concentration is calculated using the following formula:

(2)\[C_{SG,local}(i,j) = \sum_{i'=1}^{n_{x}} \sum_{j'=1}^{n_{y}} \frac{E_{SG,local}(i',j')}{U(i,j,i',j')} I\left(r(i,j,i',j'),h(i',j'),z(i,j)\right)\]

Where \(E_{SG,local}\) are the emissions attributed to the subgrid, \(U\) is the wind speed, dependent on both source and receptor subgrid values, \(n_{x}\) and \(n_{y}\) define the extent of the subgrid calculation window, \(I\) is a function determining the dispersion intensity of the emission source with horizontal spatial vectors \(r\) between the receptor grid points \((i,j)\) and the source grid points \((i',j')\) at receptor height \(z\) and source height \(h\). The contribution from each proxy emissions subgrid within the calculation window is summed at the receptor subgrid, which is centered within this window (TODO: Figure 1).

When using the independent emissions method, the local subgrid emissions, \(E_{SG,local}\), are specified directly. In contrast, when using the emission redistribution method, \(E_{SG,local}\) is derived from the EMEP emission grid, \(E_{G}(I,J)\), and the proxy data for emissions, \(P_{emission}\), normalized within the EMEP grid as follows:

(3)\[E_{SG,local}(i',j') = E_{G}(I,J) \frac{P_{emission}(i',j')}{ \sum_{i'=1}^{n_{x}} \sum_{j'=1}^{n_{y}} P_{emission}(i',j') }\]

EMEP non-local contribution calculation

(4)\[C_{G,local}(I,J,I_{lf},J_{lf},s) = LF(I,J,I_{lf},J_{lf},s) \; C_{G}(I,J)\]

(5)\[C_{G,nonlocal}(I,J) = C_{G}(I,J) - \left\{ \sum_{I'=-n_{mw}/2}^{+n_{mw}/2} \sum_{J'=-n_{mw}/2}^{+n_{mw}/2} \sum_{s=1}^{n_{source}} \; C_{G,local}(I,J,I',J',s) \right\}\]

Moving windows calculations

(6)\[C_{G,local}(i,j,s) = \sum_{I'=-n_{mw}/2}^{+n_{mw}/2} \sum_{J'=-n_{mw}/2}^{+n_{mw}/2} C_{G,local}(I,J,I',J',s) \; w(i,j,I',J',s)\]

Two methods are available for calculating these weights: area-based and emission-based. In area-based weighting, \(w_{a}\), only the area fraction of the EMEP grid that is within the moving window for the receptor subgrid is included:

(7)\[w_{a}(i,j,I',J') = \frac{\{ a(i,j) \cap A(I',J') \}}{A(I',J')}\]

Where \(A(I',J')\) is the area of the EMEP grid and \(a(i,j)\) is the area of the moving window centered at the receptor subgrid point.

Emission-based weighting requires high-resolution emission data that are compatible with the EMEP grid dimensions. In this case, \(w_{e}\) is determined by the fraction of the total subgrid emissions within the moving window and within the EMEP grid:

(8)\[w_{e}(i,j,I',J',s) = \frac{ \sum e(i,j,I',J',s) :\in \{ a(i,j) \cap A(I',J') \} } { \sum e(i,j,I',J',s) :\in \{ A(I',J') \} }\]

Where the numerator is the sum of emissions, \(e\) within the intersection of \(a(i,j)\) and \(A(I',J')\) and the demoninator is the sum of emissions within \(A(I',J')\).

(9)\[\]

(10)\[C_{SG,nonlocal}(i,j) = \frac{1}{n_{source}} \sum_{s=1}^{n_{source}} \left\{ (C_{G,local}(i,j,s) + C_{G,nonlocal}(i,j,s)) \right\} - \sum_{s=1}^{n_{source}} \left\{ C_{G,local}(i,j,s) \right\}\]

Processes and parameterization

Gaussian plume model for hourly calculations

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

(11)\[C(x,y,z) = \frac{Q}{U} I(x,y,z)\]

(12)\[I(x,y,z) = \frac{1}{2 \pi \sigma_{y} \sigma_{z}} \mathrm{exp}\left(\frac{-y^{2}}{2 \sigma_{y}^{2}}\right) \sum_{i=1}^{i=6} \left\{ \mathrm{exp}(\frac{-(z - h_{i})^{2}}{2 \sigma_{z}^{2}}) \right\}\]

(13)\[h_{i} = [ h_{emis}, -h_{emis}, 2H - h_{emis}, 2H + h_{emis}, -2H + h_{emis}, -2H - h_{emis} ]\]

(14)\[I(x,y,z) = \frac{1}{\sqrt{2 \pi} \sigma_{y} H} \mathrm{exp}\left(\frac{-y^{2}}{2 \sigma_{2}^{2}}\right)\]

(15)\[\sigma_{z}(t) = \sigma_{z0} + \sqrt{2 K_{z}(z) t f_{t}}\]

(16)\[\sigma_{y}(t) = \sigma_{y0} + \sqrt{2 K_{y}(z) t f_{t}}\]

(17)\[f_{t} = 1 + \left(\frac{\tau_{l}}{t} \mathrm{exp}\left(- \frac{t}{\tau_{l}}\right) - 1\right)\]

Where \(\tau_{l}\) is the Lagrangian integral timescale, calculated following Hanna[3] as:

(18)\[\tau_{l} = 0.6 \frac{\mathrm{max}(z_{emis}, z_{\tau min})}{u_{*}}\]

Where \(z_{emis}\) is the emission height, \(z_{\tau min}\) is the minimum emission height (2 meters), and \(u_{*}\) is the friction velocity.

(19)\[t = \frac{\mathrm{max}(x_{min},x)}{U(z)}\]

Where \(x_{min}\) is half a subgrid, \(U(z)\) is the vertical wind speed profile, and \(x\) is the down-plume distance.

(20)\[z_{cm} = \frac{\int_{0}^{H} z I(z) dz}{\int_{0}^{H} I(z) dz}\]

Solving this integral analytically yields:

(21)\[ \begin{align}\begin{aligned}z_{cm} = &\frac{\sigma_{z}}{\sqrt{2\pi}} \sum_{i=1}^{i=6} \mathrm{exp}\left( \frac{-h_{i}^{2}}{2\sigma_{z}^{2}} \right) - \mathrm{exp}\left( \frac{-(H - h_{i})^{2}}{2\sigma_{z}^{2}} \right)\\&+ \frac{h_{i}}{2} \left( \mathrm{erf}\left( \frac{h_{i}}{\sqrt{2}\sigma_{z}}\right) + \mathrm{erf}\left( \frac{(H - h_{i})}{\sqrt{2}\sigma_{z}} \right) \right)\end{aligned}\end{align} \]

For the the well-mixed case where \(\sigma_{z} > 0.9H\):

(22)\[z_{cm} = 0.5H\]

The vertical wind profile is calculated based on decreasing turbulent shear with height following Gryning et al.[4]:

(23)\[ \begin{align}\begin{aligned}U(z) = \frac{u_{* 0}}{\kappa} \left( \mathrm{log}\left( \frac{z}{z_{0}} \right) - \psi_{m} + \kappa \frac{z}{z_{l}} \left( 1 - \frac{z}{2H} \right) - \frac{z}{H} \left( 1 + b \frac{z}{2L} \right) \right) \; \mathrm{for} \; L \ge 0\\U(z) = \frac{u_{* 0}}{\kappa} \left( \mathrm{log}\left( \frac{z}{z_{0}} \right) - \psi_{m} + \kappa \frac{z}{z_{l}} \left( 1 - \frac{z}{2H} \right) - \frac{z}{H} \frac{((az - L) \phi_{m} + L)}{a(p + 1)} \right) \; \mathrm{for} \; L < 0\end{aligned}\end{align} \]

(24)\[K_{y}(z) = \frac{\sigma_{v}(z)^{2}}{\sigma_{w}(z)^{2}} K_{z}(z)\]

Rotationally symmetric Gaussian plume model for annual mean calculations

(25)\[\sigma_{y,z} = \sigma_{0(y,z)} + a_{y,z} x^{b_{y,z}}\]

(26)\[ \begin{align}\begin{aligned}I_{rot}(r,z) =& \frac{1}{\pi \sqrt{2\pi}r\varepsilon_{z}\sqrt{1 + B}} \mathrm{erf}\left( \frac{\pi \sqrt{1 + B}}{\sqrt{2}\varepsilon_{\theta}} \right)\\&\sum_{i=1}^{i=6} \left\{ \mathrm{exp} \left( \frac{-(z - h_{i})^{2}}{2\varepsilon_{z}^{2}} \right) \right\}\end{aligned}\end{align} \]

where

(27)\[\varepsilon_{z} = \sigma_{0z} + a_{z} r^{b_{z}}\]

(28)\[\varepsilon_{\theta} = \frac{1}{r} \left( \sigma_{0y} + a_{y} r^{b_{y}} \right)\]

(29)\[B = -\varepsilon_{\theta}^{2} \left( \frac{b_{z}(\varepsilon_{z} - \sigma_{0z})}{r \varepsilon_{\theta}} + \frac{b_{y}(r\varepsilon_{\theta} - \sigma_{0y})}{\varepsilon_{z}} \right)\]

(30)\[ \begin{align}\begin{aligned}I_{rot}(r,z) =& \frac{1}{2 \pi r \varepsilon_{z} \varepsilon_{\theta}} \left( 1 - \frac{\pi^{2} (1 + B)}{6\varepsilon_{\theta}^{2}} + \frac{\pi^{4}(1 + B)^{2}}{40\varepsilon_{\theta}^{4}} \right)\\&\sum_{i=1}^{i=6} \left\{ \mathrm{exp} \left( \frac{-(z - h_{i})^{2}}{2\varepsilon_{z}^{2}} \right) \right\} \; \mathrm{for} \; B < -1\end{aligned}\end{align} \]

Initial dispersion

(31)\[\sigma_{0y} = \sigma_{init,y} + 0.8 \frac{\Delta y}{2}\]

NO₂ chemistry for hourly means

(32)\[ \begin{align}\begin{aligned}\mathrm{NO + O_{3} \to NO_{2} + O_{2}} \qquad \mathrm{(a)}\\\mathrm{NO_{2}} + hv \to \mathrm{NO + O} \qquad \mathrm{(b)}\\\mathrm{O_{2} + O} + M \to \mathrm{O_{3}} + M \qquad \mathrm{(c)}\end{aligned}\end{align} \]

(33)\[\frac{\mathrm{d[NO_{2}]}}{\mathrm{d}t} = k_{1} \mathrm{[NO][O_{3}]} - J \mathrm{[NO_{2}]}\]

(34)\[k_{1} = 1.4 \times 10^{-12} \mathrm{exp} \left( \frac{-1310}{T_{air}} \right) (\mathrm{cm^{3} \; s^{-1}})\]

References

[1]

Bruce Rolstad Denby, Michael Gauss, Peter Wind, Qing Mu, Eivind Grøtting Wærsted, Hilde Fagerli, Alvaro Valdebenito, and Heiko Klein. Description of the uEMEP_v5 downscaling approach for the EMEP MSC-W chemistry transport model. Geoscientific Model Development, 13(12):6303–6323, December 2020. doi:10.5194/gmd-13-6303-2020.

[2]

P. Wind, B. Rolstad Denby, and M. Gauss. Local fractions – a method for the calculation of local source contributions to air pollution, illustrated by examples using the emep msc-w model (rv4_33). Geoscientific Model Development, 13(3):1623–1634, 2020. doi:10.5194/gmd-13-1623-2020.

[3]

Steven R. Hanna. Lagrangian and Eulerian Time-Scale Relations in the Daytime Boundary Layer. Journal of Applied Meteorology, 20(3):242–249, March 1981. doi:10.1175/1520-0450(1981)020<0242:LAETSR>2.0.CO;2.

[4]

Sven-Erik Gryning, Ekaterina Batchvarova, Burghard Brümmer, Hans Jørgensen, and Søren Larsen. On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer. Boundary-Layer Meteorology, 124(2):251–268, July 2007. doi:10.1007/s10546-007-9166-9.