A quasi-experimental coastal region eddy diffusivity applied in the APUGRID model

In this study, Taylor statistical diffusion theory and sonic anemometer measurements collected at 11 levels on a 140 m high tower located in a coastal region in southeastern Brazil have been employed to obtain quasi-empirical convective eddy diffusivity parameterizations in a planetary boundary layer (PBL). The derived algebraic formulations expressing the eddy diffusivities were introduced into an Eulerian dispersion model and validated with Copenhagen tracer experiments. The employed Eulerian model is based on the numerical solution of the diffusion–advection equation by the fractional step/locally one-dimensional (LOD) methods. Moreover, the semi-Lagrangian cubic-spline technique and Crank–Nicolson implicit scheme are considered to solve the advection and diffusive terms. The numerical simulation results indicate that the new approach, based on these quasi-experimental eddy diffusivities, is able to reproduce the Copenhagen concentration data. Therefore, the new turbulent dispersion parameterization can be applied in air pollution models.

1 Introduction

Eulerian models are powerful tools to study and investigate the air pollution dispersion in the planetary boundary layer (PBL) (Hanna et al., 1982; Tirabassi, 2009; Zannetti, 2013). These models are based on the solution of the classical advection–diffusion equation, containing the turbulent eddy diffusivities, which provide the simulated contaminant concentration data (Batchelor, 1949; Pasquill and Smith, 1983).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial t}}& +\stackrel{\mathrm{‾}}{u}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial x}}+\stackrel{\mathrm{‾}}{v}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial y}}+\stackrel{\mathrm{‾}}{w}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial z}}\\ & =-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u^{\prime }{C}^{\prime }}}{\partial x}}-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v^{\prime }{C}^{\prime }}}{\partial y}}-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{w^{\prime }{C}^{\prime }}}{\partial z}}+S,\end{array}\end{array}$$

where C is the concentration and S is a source term. These Eulerian models describe the concentration turbulent fluxes as the gradient of the mean concentration employing the eddy diffusivities (K theory):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{u^{\prime }{C}^{\prime }}& =-K_x{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial x}},\\ \stackrel{\mathrm{‾}}{v^{\prime }{C}^{\prime }}& =-K_y{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial y}},\\ \stackrel{\mathrm{‾}}{w^{\prime }{C}^{\prime }}& =-K_z{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial z}},\end{array}\end{array}$$

where K_x, K_y, and K_z are the eddy diffusivities in the x, y, and z directions and $\stackrel{\mathrm{‾}}{u}$, $\stackrel{\mathrm{‾}}{v}$, and $\stackrel{\mathrm{‾}}{w}$ represent the longitudinal, lateral, and vertical mean wind components, respectively. Thus, Eq. (1) can be written in the form

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial t}}& +\stackrel{\mathrm{‾}}{u}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial x}}+\stackrel{\mathrm{‾}}{v}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial y}}+\stackrel{\mathrm{‾}}{w}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial z}}={\displaystyle \frac{\partial }{\partial x}}\left(K_x{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial x}}\right)\\ & +{\displaystyle \frac{\partial }{\partial y}}\left(K_y{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(K_z{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial z}}\right)+S.\end{array}\end{array}$$

From the numerical point of view, to solve Eq. (3) it is required to provide the wind and turbulence physical description. For the turbulent diffusion we need to specify K_x, K_y and K_z. These turbulent parameters with dimensions of length times velocity therefore describe the eddy size and eddy velocity (Panofsky and Dutton, 1984).

Most of the eddy diffusivities employed in current operational dispersion models are based on PBL similarity theories (Leelőssy et al., 2014). However, a better description of the turbulent properties, associated with eddy diffusivities, is based on direct measurements of wind data with high vertical resolution (Martins et al., 2018).

The coastal internal boundary layers (CIBLs) are generated by differences in surface temperature and aerodynamic roughness occurring between land and water atmospheric environments. Considering that a large number of power plants and industrial complexes and hence polluting installations are constructed in coastal regions, it is necessary to obtain CIBL turbulent parameters that are employed in dispersion models to describe the coastal air pollution. The growing interest in the dispersion issues regarding pollutant emission in coastal areas requires knowledge of the turbulent structure of the planetary boundary layer in this region. However, the characteristics of the turbulence in these boundary layers vary complexly in space and time due to the sudden changes in the surface characteristics, as heat flux and roughness, in the sea–land interface. In the occurrence of sea breeze, the stably stratified air mass over the water reaches the coast and starts to be heated by the land surface. Thus, a convective boundary rises from the surface, developing a thermal internal boundary layer (TIBL) that increases in height as it advances over the land. The TIBL is topped by a stably stratified inversion layer that affects the atmospheric diffusion in coastal regions. Therefore, to improve the response of the dispersion models, it is necessary to provide a truthful description of the turbulence through the TIBL. In this sense, several observational experiments are performed using airborne, tethered balloons and fixed mast measurement techniques (Smedman and Hoegstroem, 1983; Ogawa and Ohara, 1985; Durand et al., 1989; Shao et al., 1991). Wind-tunnel experiments and numerical simulations are found in Hara et al. (2009).

In this present study, we use eddy diffusivities that were derived from the observations of the turbulent wind components (u, v, w) in a convective CIBL to simulate the dispersion of contaminants released from an elevated continuous point source in a coastal region. The turbulent observations were performed at a 140 m micrometeorological tower positioned 240 m north of a natural gas power plant and 4 km southwest of the ocean coastal environment in the city of Linhares (southeastern Brazil). The turbulent wind data were obtained from high-frequency measurements (10 Hz) accomplished by tridimensional sonic anemometers at heights of 1, 2, 5, 9, 20, 37, 56, 75, 94, 113 and 132 m (Martins et al., 2018). Therefore, the study of Martins et al. (2018) employs these measurements, the turbulent energy spectra and a few mathematical relations to determine turbulent dispersion parameters (Taylor statistical diffusion theory, Degrazia et al., 2000, 2001).

Differently, of previous studies in which the vertical profiles of turbulent parameters have been calculated using surface observations to throughout the similarity-based relationship, our eddy diffusivities were locally calculated from the detailed measurements accomplished along the entire vertical extension occupied by the surface internal boundary layer. As a consequence, they can be called quasi-experimental eddy diffusivities. The aim of this work is to obtain algebraic formulations from the fitting curves that reproduce the observed vertical profile of these quasi-experimental eddy diffusivities. As a test and to evaluate the quasi-experimental eddy diffusivities for a convective CIBL, we substitute these turbulent diffusion parameters into Eq. (3) to simulate the contaminant concentration originating from an elevated continuous point source in a coastal environment. The simulated concentrations are compared to those measured in the Copenhagen diffusion experiments.

From the point of view of originality and novelty, the present development, from some asymptotic equations and detailed turbulent spectral observations of the surface coastal internal boundary layer, provides a general methodology for obtaining algebraic expressions that reliably represent the eddy diffusivities in the coastal internal boundary layer.

2 Eulerian grid-dispersion model

In this section, to simulate the contaminant concentrations using Eq. (3), we present the Eulerian grid-dispersion model proposed by Rizza et al. (2003), the so-called APUGRID. The APUGRID model employs fractional step/locally one-dimensional (LOD) methods (Yanenko, 1971; McRae et al., 1982; Mařcuk, 1984) to solve the diffusion–advection equation. In the LOD numerical method, Eq. (3) is separated into time-dependent equations, each one locally one-dimensional (LOD) (Yanenko, 1971; Rizza et al., 2003; Yordanov et al., 2006). As a consequence,

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{C}}{\partial t}}={\mathrm{\Lambda }}_x\stackrel{\mathrm{‾}}{C}+{\mathrm{\Lambda }}_y\stackrel{\mathrm{‾}}{C}+{\mathrm{\Lambda }}_z\stackrel{\mathrm{‾}}{C},\end{array}$$

where ${\mathrm{\Lambda }}_x=-\stackrel{\mathrm{‾}}{u}\frac{\partial }{\partial x}+\frac{\partial }{\partial x}\left(K_x\frac{\partial }{\partial x}\right)$, ${\mathrm{\Lambda }}_y=-\stackrel{\mathrm{‾}}{v}\frac{\partial }{\partial y}+\frac{\partial }{\partial y}\left(K_y\frac{\partial }{\partial y}\right)$ and ${\mathrm{\Lambda }}_z=-\stackrel{\mathrm{‾}}{w}\frac{\partial }{\partial z}+\frac{\partial }{\partial z}\left(K_z\frac{\partial }{\partial z}\right)$.

Employing Crank–Nicolson time integration (McRae et al., 1982; Yordanov et al., 2006; Rizza et al., 2010), we obtain that

$$\begin{array}{}\text{(5)}& {\stackrel{\mathrm{‾}}{C}}^{n+\mathrm{1}}=\prod _{j=\mathrm{1}}^{\mathrm{3}}{\left[\mathbf{I}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}\right]}^{-\mathrm{1}}{\left[\mathbf{I}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\mathrm{\Lambda }}_j\right]}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{C}}^n,\end{array}$$

with I being the unity matrix and Δt the time step. The second-order accuracy can be obtained following Rizza et al. (2010) by

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{‾}}{C}}^n=\prod _{i=\mathrm{1}}^{\mathrm{3}}\left[A_{i}F{D}_i\right]{\stackrel{\mathrm{‾}}{C}}^{n+\mathrm{1}}{\stackrel{\mathrm{‾}}{C}}^{n+\mathrm{1}}=\prod _{i=\mathrm{1}}^{\mathrm{3}}\left[D_i{A}_{i}F\right]{\stackrel{\mathrm{‾}}{C}}^n.\end{array}$$

In Eq. (6) $A_i=\stackrel{\mathrm{‾}}{u_i}\frac{\partial }{\partial i}$, $D_i=\frac{\partial }{\partial i}\left(K_i\frac{\partial }{\partial i}\right)$, with $i=x,y,z$ and F representing the filter operation.

The advection terms were solved by employing a quasi-Lagrangian cubic-spline technique (Long and Pepper, 1981), and the numerical model stability is carried out by the Courant–Friedrichs–Lewy condition (cfl):

$$\begin{array}{}\text{(7)}& \mathrm{cfl}=\stackrel{\mathrm{‾}}{U}{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}},\end{array}$$

$\stackrel{\mathrm{‾}}{U}$ being the mean wind speed and Δx the grid spacing with the stability condition cfl⩽1 satisfied.

In order to calculate the concentration advective transport by the mean wind speed, we use the wind speed profile described by the following similarity law (Berkowicz et al., 1986):

$$\begin{array}{}\text{(8)}& \stackrel{\mathrm{‾}}{U}\left(z\right)={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\left[\log (z/z_{\mathrm{0}})-{\mathit{\psi }}_m(z/L)+{\mathit{\psi }}_m(z_{\mathrm{0}}/L)\right].\end{array}$$

Equation (8) is valid for z<z_b, where z_b=0.1z_i, where z_i is the convective boundary layer height, u_* is the friction velocity, κ is the von Kármán constant and z₀ is the surface roughness. ψ_m is a stability function defined as

$$\begin{array}{}\text{(9)}& {\mathit{\psi }}_m=\mathrm{2}\ln \left[{\displaystyle \frac{\mathrm{1}+A}{\mathrm{2}}}\right]+\ln \left[{\displaystyle \frac{\mathrm{1}+A^{\mathrm{2}}}{\mathrm{2}}}\right]-\mathrm{2}{\tan }^{-\mathrm{1}}\left(A\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},\end{array}$$

with $A={\left(\mathrm{1}-\mathrm{16}z/L\right)}^{\mathrm{1}/\mathrm{4}}$ and L the Obukhov length. For z>z_b the wind profile is the wind speed at z=z_b.

Quasi-empirical eddy diffusivity models: evaluation in APUGRID

The eddy diffusivities can be found by the following relationship:

$$\begin{array}{}\text{(10)}& K_i={\mathit{\sigma }}_i^{\mathrm{2}}{T}_{\mathrm{L}i},\end{array}$$

where ${\mathit{\sigma }}_i^{\mathrm{2}}$ is the turbulent velocity variance quantifying the turbulence mixing degree and T_Li is the decorrelation local timescale that takes into account the characteristic time in which a fluid control volume maintains its motion in a particular direction (Hinze, 1975).

To obtain the Lagrangian K_i from the Eulerian measurements, the relation between the Lagrangian T_Li and Eulerian decorrelation T_Ei timescales are employed (Hanna, 1981; Degrazia and Anfossi, 1998):

$$\begin{array}{}\text{(11)}& {\mathit{\beta }}_i={\displaystyle \frac{T_{\mathrm{L}i}}{T_{\mathrm{E}i}}}={\displaystyle \frac{\mathrm{0.55}\stackrel{\mathrm{‾}}{U}}{{\mathit{\sigma }}_i}},\end{array}$$

where

$$\begin{array}{}\text{(12)}& T_{\mathrm{E}i}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_i^{\mathrm{2}}}}\left[\stackrel{\mathrm{‾}}{u_i^{\prime }+u_i^{\prime }(t+\mathit{\tau })}\right]\mathrm{d}\mathit{\tau },\end{array}$$

where $u_i^{\prime }$ is the wind speed turbulent fluctuation and τ is the temporal lag. Recently, Martins et al. (2018) used Eqs. (10), (11) and (12) to derive experimental vertical profiles for K_x, K_y and K_z. To obtain such profiles, 1 h observation wind velocity time series intervals are tested for quality control requirements. Unstable conditions were considered to be daytime time series with $-\mathrm{150}\le L<\mathrm{0}$. From a total of 4 months of observations (August–November 2016), 343 1 h unstable intervals are retained. The variances and timescale profiles used to estimate the K_i vertical profiles are obtained by averaging all 343 individual profiles.

From this set of eddy diffusivity vertical profiles, we use the best fit curves approach to obtain the following simple algebraic formulations:

$$\begin{array}{}\text{(13)}& {\displaystyle }{K}_x\left(z\right)=\mathrm{7.83}\times {\mathrm{10}}^{\mathrm{1}}-\mathrm{6.42}\times {\mathrm{10}}^{\mathrm{1}}z,\text{(14)}& {\displaystyle }{K}_y\left(z\right)=\mathrm{7.35}\times {\mathrm{10}}^{\mathrm{1}}\log \left(z\right)+\mathrm{4.25}z-\mathrm{3.73}\times {\mathrm{10}}^{-\mathrm{2}}{z}^{\mathrm{2}},\text{(15)}& {\displaystyle }{K}_z\left(z\right)=\mathrm{8.34}\times {\mathrm{10}}^{-\mathrm{1}}z.\end{array}$$

In order to test these eddy diffusivities, we perform contaminant concentration simulations and compare the simulated data with the Copenhagen tracer dispersion experiments (Gryning and Lyck, 1984). The tracer sulfur hexafluoride (SF6) used in the Copenhagen dispersion experiments was released at a height of 115 m from the TV tower in Gladsaxe (Copenhagen) and the ground-level contaminant concentrations were measured at 3 arcsec located at a distance of 2000 to 6000 m from the elevated continuous point source. The experiment site is limited by the Øresund coast, approximately 7 km east of the TV tower. Therefore, the turbulent effects acting on the tracer dispersion are characteristic of the CIBL. The width of Øresund, the water portion separating Denmark and Sweden, is about 20 km. On the western side of Øresund lies Copenhagen with its urban area. This area has high surface roughness due to the urban character. Thusly, a turbulent environment occurs in a region with relatively cold water and a warm land surface. As a consequence, the turbulent structure acting on the tracer dispersion can be considered to be one present in the coastal inner boundary layer. Meteorological parameters for the Copenhagen runs are shown in Table A1 in Appendix A, ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{115}\mathrm{m}}$ and ${\stackrel{\mathrm{‾}}{U}}_{\mathrm{10}\mathrm{m}}$ being the mean wind velocity measured at 115 and 10 m, respectively, σ_w the vertical wind velocity variance and z_i the convective boundary layer depth. Although some Copenhagen dispersion experiments occurred in quasi-neutral conditions, the L parameter was negative. The presence of a slightly convective stratified boundary layer can be seen in u and v turbulent energy spectra (Kaimal et al., 1972; Martins et al., 2018). In this situation, a structure that contains two peaks can be observed in spectral curves: one low-frequency peak and one high-frequency peak. This reflects the impact of the larger convective eddies on the turbulent structure (Garratt, 1992).

The choice of the Copenhagen experiment was motivated by the fact that the region in which the experiment occurred is located near the Øresund coast. The eddy diffusivities obtained from the Linhares ocean coastal environment are empirical. Therefore, it is reasonably relevant that they can be employed to simulate concentration data in another different coastal area such as the suburbs of Copenhagen. In this aspect, it can be said that although the Copenhagen data set is composed of a limited number of runs, this comparison is suitable only for preliminary validation.

Therefore, we expect that our eddy diffusivities obtained at a coastal site localized in southeastern Brazil are adequate for reproducing contaminant data in coastal regions.

In Table A2, the predicted crosswind-integrated concentrations obtained from the APUGRID model are compared, with Copenhagen diffusion experiments, over the different distances of the release point source. The performance of the APUGRID model employing the quasi-experimental eddy diffusivities as given by Eqs. (13), (14) and (15) to simulate the Copenhagen observation data can also be evaluated by analyzing the results shown in Table A3 and Fig. A1 in Appendix A.

Table A3 exhibits Hanna's statistical indices, which are commonly used to calibrate air pollution dispersion models. Such indices are defined as

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathrm{\text{normalized mean square error (NMSE)}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\left(C_{\mathrm{o}}-C_{\mathrm{p}}\right)}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{C_{\mathrm{o}}{C}_{\mathrm{p}}}}},\text{(17)}& {\displaystyle }\text{correlation coefficient}\left(R\right)={\displaystyle \frac{\left(C_{\mathrm{o}}-\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}\right)\left(C_{\mathrm{p}}-\stackrel{\mathrm{‾}}{C_{\mathrm{p}}}\right)}{{\mathit{\sigma }}_{\mathrm{o}}{\mathit{\sigma }}_{\mathrm{p}}}},\text{(18)}& {\displaystyle }\text{fractional bias (FB)}={\displaystyle \frac{\left(\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}-\stackrel{\mathrm{‾}}{C_{\mathrm{p}}}\right)}{\mathrm{0.5}\left(\stackrel{\mathrm{‾}}{C_{\mathrm{o}}}+\stackrel{\mathrm{‾}}{C_{\mathrm{p}}}\right)}},\text{(19)}& {\displaystyle }\text{fractional standard deviations (FS)}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{o}}-{\mathit{\sigma }}_{\mathrm{p}}}{\mathrm{0.5}\left({\mathit{\sigma }}_{\mathrm{o}}{\mathit{\sigma }}_{\mathrm{p}}\right)}},\end{array}$$

where C_p is the predicted concentration, C_o is the observed concentration, σ_p is the predicted standard deviation, σ_o is the observed standard deviation, and the overbar represents an averaged value.

The observed and predicted scatter diagram of concentrations in Fig. A1 demonstrates that the simulated concentration reproduces fairly well the measured concentration data. Furthermore, the statistical analysis of the results (Table A3) shows good agreement between the results of the proposed approach with the experimental ones. The indices are found within an acceptable interval, with NMSE (normalized mean square error), FB (fractional bias) and FS (fractional standard deviation) close to zero and R (correlation coefficient) near to one. These statistical indices show that the empirical eddy diffusivities obtained at a Brazilian coastal site can be used to simulate contaminant dispersion in other coastal areas. Thus, the present development based on an analysis of high-resolution turbulence data from an elevated micrometeorological tower provides suitable eddy diffusivities that describe the turbulent transport patterns in a CIBL.

3 Conclusions

The Eulerian operational air dispersion models that simulate contaminant observed concentration data need to incorporate into their formulation the characteristics of the PBL turbulent diffusion process. To accomplish this parameterization they use turbulent transport terms known as eddy diffusivities. These turbulent parameters represent approximate quantities which intend to reproduce the complex natural dispersive effects. In this study, algebraic expressions for quasi-experimental convective eddy diffusivities for a coastal site are derived. The derivation employed the Taylor statistical diffusion theory and sonic anemometer observations with high vertical resolution in a CIBL.

The complexity of the subject does not allow a direct confrontation between experiment and model. However, utilizing the APUGRID Eulerian dispersion model and a concentration data set of dispersion experiments performed in a CIBL, the derived eddy diffusivities have been tested and validated. The comparison results show that there is fairly good agreement between simulated and measured concentrations. As a consequence, the results provided in this investigation are encouraging. Thus, the new eddy diffusivities for a coastal site may be suitable for applications in regulatory air pollution modeling.

Appendix A

Figure A1 Scatter diagram between Copenhagen observed (C_o∕Q) and predicted (C_p∕Q) ground-level crosswind-integrated concentrations normalized by the emission rate.

Table A1 Meteorological conditions during the Copenhagen dispersion experiments.

Table A2 Observed C_o and predicted C_p crosswind-integrated concentrations normalized by the emission rate (Q) for Copenhagen experiments.

Table A3 Statistical evaluation of the APUGRID model employing the quasi-experimental eddy diffusivities.