Intercomparison of two K-schemes: Local versus non-local in calculating concentrations of pollutants in chemical and air-quality models

Environmental Modelling & Software 22 (2007) 1685e1689 www.elsevier.com/locate/envsoft

Short communication

Intercomparison of two K-schemes: Local versus non-local in calculating concentrations of pollutants in chemical and air-quality models D.T. Mihailovic a,*, K. Alapaty b,1 a

Faculty of Agriculture, University of Novi Sad, Trg Dositeja Obradovica 8, 21000 Novi Sad, Serbia and Montenegro b Division of Atmospheric Sciences, National Science Foundation, Arlington, VA 22230, USA Received 10 October 2006; received in revised form 9 March 2007; accepted 9 March 2007 Available online 23 April 2007

Abstract In this paper, we studied the performance of a local and non-local scheme for vertical diffusion in the atmospheric boundary layer and their impact on the concentration of pollutants calculated with chemical and air-quality models. In the local diffusion scheme, the vertical eddy-diffusivity is determined independently at each point based on local vertical gradients of wind and potential temperature. The non-local scheme determines an eddy-diffusivity profile based on a diagnosed boundary layer height and a turbulent vertical scale. Depending on the boundary layer stability regime, a modification was introduced into the non-local scheme to facilitate different values of the parameter in the equation for eddy-diffusivity. To examine the performance of the schemes, simulated and measured concentrations of a pollutant (NO2), believed to be one of the most affected ones by the processes in the atmospheric boundary layer, were compared for the years 1999, 2001 and 2002. The comparison was made for the whole domain used in simulations performed by the chemical EMEP Unified model (version UNI-ACID, rv2.0) where schemes were incorporated. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Atmospheric boundary layer; Non-local closure model; Turbulent kinetic energy; Atmospheric chemistry; Environmental modelling

1. Introduction One of the well-known issues regarding local-closure atmospheric boundary layer (ABL) schemes is their inability to produce well-mixed layers in the ABL during convective conditions. Holtslag and Boville (1993), using the NCAR Community Climate Model (CCM2), studied a classic example of artifacts resulting from the deficiencies in the first-order closure schemes. To alleviate problems associated with the general first-order eddy-diffusivity K-schemes, they proposed a non-local K-scheme. Hong and Pan (1996) presented an enhanced version of the Holtslag and Boville (1993) scheme. In * Corresponding author. Tel./fax: þ381 21 6350 552. E-mail address: [email protected] (D.T. Mihailovic). 1 Present address: Office of Biological and Environmental Research, Office of Science, U.S. Department of Energy, USA. 1364-8152/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2007.03.002

this scheme the friction velocity scale (u*) is used as a closure in their formulation. However, for moderate to strong convective conditions, u* is not a representative scale (Alapaty and Alapaty, 2001). Rather, the convective velocity (w*) scale is suitable as used by Hass et al. (1991) in simulation of a wet deposition case in Europe by the European Acid Deposition Model (EURAD). Depending on the magnitude of the scaling parameter h/L (h is depth of the ABL, and L is MonineObukhov length), either u* or w* is used in many other formulations. Notice that this approach may not guarantee continuity between the alternate usage of u* and w* in estimating K-eddy diffusivity. Also, in most of the local-closure schemes the coefficient of vertical eddy-diffusivity for moisture is assumed to be equal to that for heat. Sometimes this assumption leads to vertical gradients in the simulated moisture fields, even during moderate to strong convective conditions in the ABL. Also, the non-local scheme considers the horizontal

D.T. Mihailovic, K. Alapaty / Environmental Modelling & Software 22 (2007) 1685e1689

1686

advection of turbulence that may be important over heterogeneous landscapes (Alapaty and Alapaty, 2001; Mihailovic et al., 2005). In this paper, we propose a modified version of a first-order non-local K-scheme (Alapaty, 2003), addressing the deficiencies discussed above. The main feature of the scheme that has been used for tests in this study is shortly described in Section 2. The results of comparison with observations obtained by the Unified EMEP model, using local and non-local diffusivity schemes for the years 1999, 2001 and 2002, are presented in Section 3.

The formulation of eddy-diffusivity by Eq. (1) depends on h. We follow Troen and Mahrt (1986) for determination of h using,

 Ric uðhÞ2 þvðhÞ2 h¼ g ; ð5Þ fqv ðhÞ  qs g T0 where Ric is a critical bulk Richardson number for the ABL, u(h) and v(h) are the horizontal velocity components at h, g=q0 is the buoyancy parameter, qs is the appropriate virtual potential temperature, and qv ðhÞ is the virtual potential temperature of air near the surface at h, respectively. For unstable conditions (L < 0), qs is given by Troen and Mahrt (1986). qs ¼ qs ðz1 Þ þ C0

2. Methodology

The starting point of approach is to consider the general form of the vertical eddy-diffusivity equation. For momentum, this equation can be written as:

 2=3  2=3 1 LE Fm ; 0:4w3 þ u3 ðh  zÞ 2 h kz

ð1Þ

ð3Þ

Following works of Zhang et al. (1996) and Moeng and Sullivan (1994) about LES (Large Eddy Simulation), Alapaty (2003) suggested how to estimate the vertically integrated mean turbulent velocity scale e that in the ABL can be written as: 1 e ¼ h

Z h pffiffiffiffiffiffiffiffi eðzÞJðzÞdz;

ð7Þ

w ¼ ððg=q0 Þwqv0 hÞ1=3 :

ð8Þ

In Eq. (6), qv ðz1 Þ is the virtual temperature at the first model level. The second term on the right-hand side of Eq. (6) represents a temperature excess, which is a measure in the lower part of the ABL. For stable conditions we use qs ¼ qv ðz1 Þ with z1 ¼ 2 m. On the basis of Eq. (5) the boundary layer height can be calculated by iteration for all stability conditions, when the surface fluxes and profiles of qv, u and v are known. The computation starts with calculating the bulk Richardson number Ri between the level qs and subsequent higher levels of the model. Once Ri exceeds the critical value, the value of h is derived with linear interpolation between the level with Ri > Ric and the level underneath. We use a minimum of 100 m for h. In Eq. (5), Ric is the value of the critical bulk Richardson number used to be 0.25 in this study. In the free atmosphere, turbulent mixing is parameterized using the formulation suggested by Blackadar (1979) in which vertical eddy diffusivities are functions of the Richardson number and wind shear in the vertical. This formulation can be written as: Km ¼ K0 þ SðklÞ2

Rc  Ri ; Rc

ð9Þ

where K0 is the background value (1 m2 s1), S is the vertical wind shear, l is the characteristic turbulent length scale (100 m), Rc is the critical Richardson number, and Ri is the Richardson number defined as:

ð2Þ

where LE characterizes the integral length scale of the dissipation rate. Here, Fm ¼ ð1  15z=LÞ1=4 is an empirical function for the unstable atmospheric surface layer (Businger et al., 1971), which is applied to both the surface and mixed layer. We used LE ¼ 2.6h which is in the range 2.53.0h suggested by Moeng and Sullivan (1994). For the stable atmospheric boundary layer we modelled the TKE profile using an empirical function proposed by Lenschow et al. (1988), based on aircraft observations,  eðzÞ z 1:75 ¼6 1 : 2 u h

1=3 ws ¼ u3 þ c1 w3 ; and

where Km is the vertical eddy-diffusivity, e is the mean turbulent velocity scale in the ABL to be determined (closure problem), k is the von Karman constant (k ¼ 0.41), z is the vertical coordinate, p is the profile shape exponent coming from the similarity theory (Troen and Mahrt, 1986; usually taken as 2), and Fm is the nondimensional function of momentum. According to Zhang et al. (1996), we use the square root of the vertically averaged turbulent kinetic energy (TKE) in the ABL as a velocity scale, in place of the mean wind speed, the closure to Eq. (1). Instead of using a prognostic approach to determine TKE, we make use of a diagnostic method. It is then logical to consider the diagnostic TKE to be a function of both u* and w*. Thus, the square root of diagnosed TKE near the surface serves as a closure to this problem (Alapaty and Alapaty, 2001). However, it is more suitable to estimate e from the profile of the TKE through the whole ABL. According to Moeng and Sullivan (1994), a linear combination of the turbulent kinetic energy dissipation rates associated with shear and buoyancy can adequately approximate the vertical distribution of the turbulent kinetic energy (TKE), e(z), in a variety of boundary layers ranging from near neutral to free convection conditions. Following Zhang et al. (1996) the TKE profile can be expressed as: eðzÞ ¼

ð6Þ

where C0 ¼ 8.5 (Holtslag et al., 1990), ws is the velocity, while wqv0 is the kinematics surface heat flux. The velocity ws is parameterized as:

2.1. Description of the TKE vertical diffusion scheme

p e kzð1  hz ; Km ¼ Fm

wqv0 ; ws

ð4Þ

0

where JðzÞ is the vertical profile function for turbulent kinetic energy as obtained by Zhang et al. (1996) based on LES studies, later modified by Alapaty (personal communication), and dz is layer thickness.

g vqv : qv S2 vz The critical Richardson number in Eq. (9) is determined as:

Ri ¼

ð10Þ

Rc ¼ 0:257ðDzÞ0:175 ;

ð11Þ

where Dz is the layer thickness (Zhang and Anthes, 1982).

2.2. Formulation of vertical diffusion in the Unified EMEP chemical model The vertical sub grid turbulent transport in the EMEP Unified model is modelled as a diffusivity effect. The local eddy-diffusivity scheme is designed following O’Brien (1970). (In further text this scheme will be refereed to the OLD one.) In the unstable case, Km is determined as:  2 hz Km ðzÞ ¼ Km ðhÞ þ ½Km ðhs Þ  Km ðhÞ þ ðz  hs Þ h  hs   v Km ðhs Þ  Km ðhÞ ðKm ðhs ÞÞ þ 2  hs  z < h; vz h  hs

ð12Þ

where hs is the height of the surface boundary layer. In the model calculation hs is set to 4% of the mixing height.

D.T. Mihailovic, K. Alapaty / Environmental Modelling & Software 22 (2007) 1685e1689

3. Tests and results 3.1. Short model description and experimental set up To compare performances of the proposed non-local TKE scheme (Eqs. (1)e(4)) and local OLD scheme (Eq. (12)), both based on the vertical eddy diffusivity formulation, in reproducing the vertical transport of pollutants in the ABL, a test was performed with the Unified EMEP chemical model (UNIT-ACID, rv2_0_9). The basic physical formulation of the EMEP model is unchanged from that of Berge and Jacobsen (1998). A polar-stereographic projection, true at 60 N and with the grid size of 50  50 km2 was used. The model domain used in simulation had (101, 91) points covering the area of whole Europe and North Africa. The s terrain-following coordinate was used with 20 levels in the vertical e from the surface to 100 hPa and with the lowest level located nearly at 92 m. The horizontal grid of the model is the Arakawa C grid. All other details can be found in Simpson et al. (2003). The Unified EMEP model uses three-hourly resolution meteorological data from the dedicated version of the HIRLAM (HIgh Resolution Limited Area Model) numerical weather prediction model with a parallel architecture (Bjorge and Skalin, 1995). The horizontal and vertical wind components are given on a staggered grid. All other variables are given in the centre of the grid. Linear interpolation between the three-hourly values is used to calculate values of the meteorological input data at each advection step. The time step used in the simulation was 600 s. In the EMEP Unified model the OLD scheme remarkably improved the vertical mixing in the ABL, particularly under stable conditions and conditions approaching free convection, compared with the scheme previously used in the EMEP Unified model. The improvement was particularly pronounced for NO2 (Fagerli and Eliassen, 2002). However, with reducing the horizontal grid size and increasing the heterogeneity of the underlying surface in the EMEP Unified model, there is a need for eddy-diffusivity scheme having a higher level of sophistication in the simulation of turbulence in the ABL. It seems that the non-local eddy-diffusivity schemes have good performance for that. Recently, Zhang et al. (2001) demonstrated some advantages of non-local over local eddy-diffusivity schemes. 3.2. Comparison with the observations The comparison of the proposed TKE and OLD schemes has been performed using simulated and measured concentrations of the pollutant NO2 since it is one of the most affected ones by the processes in the ABL layer. The simulations were done for the years 1999, 2001 and 2002 in the months when the convective processes are dominant in the ABL (AprileSeptember). The station recording NO2 in air (mg (N) m3) concentration was considered for comparison when measurements were available for at least 75% of days in a year [1999 (80 stations), 2001 (78) and 2002 (82)]. We have calculated the bias on the monthly basis as (M  O)/O  100% where M and O denote the modelled and observed values, respectively. The comparison of

1687

the modelled and observed NO2 in air (mg (N) m3) concentrations and corresponding biases for both schemes are shown in Fig. 1. The values used in calculations were averaged over the whole domain of integration. It can be seen that both schemes underestimate the observations. However, for all considered months, NO2 concentrations calculated with the TKE scheme are in general higher and closer to the observations than those obtained by the OLD scheme (of the order of 10%). Correspondingly, the bias of the TKE scheme is lower than the OLD scheme. To quantify the simulated values of the both schemes we have performed an error analysis of the NO2 concentration outputs NO2 based on a method discussed in Pielke (2002). Following that study, we computed several statistical quantities as follows: P P ^ 2 b i Þ2 =N1=2 ; nBR ¼ f N ½ðGi GÞðG b i GÞ n ¼ ½ Ni¼1 ðGi G i¼1 P P ^ 2 b i GÞ =Ng1=2 ; h ¼ ½ Ni¼1 ðGi GÞ2 =N1=2 , and b h ¼ ½ Ni¼1 ðG =N1=2 . Here, G is the variable of interest (aforementioned variables in this study) while N is the total number of data. An overbar indicates the arithmetic average, while a caret refers to an observation. The absence of a caret indicates a simulated value; n is the rmse, while nBR is rmse after a bias is removed. Rootmean-square errors give a good overview of a data set, with large errors weighted more than many small errors. The standard deviations in the simulations and the observations are given by h and b h. A rmse that is less than the standard deviation of the observed value indicates skill in the simulation. Moreover, the values of h and b h should be close if the prediction is to be considered realistic. The statistics gave the following values: TKE (n¼0.548, nBR ¼ 0.293, h ¼ 0.211, b h ¼ 0.147) and OLD (n ¼ 0.802, nBR ¼ 0.433, h ¼ 0.303, b h ¼ 0.147). A comparison of h and b h shows that difference between them is evidently smaller with the TKE scheme. We realise that the results depend on the different interactions of the two diffusion schemes with the other part of the model. Since no attempt was made to calibrate the diffusion schemes to give the best possible simulations, the present comparisons can only serve as illustrations of the impact.

3.3. Conclusions and recommendation In this study we have considered the impact of a local and a non-local scheme for vertical diffusion in the atmospheric boundary layer calculated with the EMEP Unified model (version UNI-ACID, rv2.0). The outputs of the local and non-local diffusion schemes (obtained from numerical experiments using 3D model with both schemes for the years 1999, 2001, 2002) have been analysed and compared for concentrations (mg m3) of a pollutant most affected by the processes in the ABL layer (NO2). The OLD scheme is used as a reference when assessing how well the TKE scheme simulates the mixing processes in the ABL. Simulated concentrations of NO2 by the TKE scheme are closer to the observations when compared to those obtained from using the OLD scheme. Further, we found that the TKE scheme is computationally more economical than the OLD scheme.

D.T. Mihailovic, K. Alapaty / Environmental Modelling & Software 22 (2007) 1685e1689 2.0

0

Year: 1999

Year: 1999 -10

1.8 1.5

Bias (M-O)/O (%)

Concentration of NO2(µg(N)m-3)

1688

OBS.

1.3

TKE

1.0 OLD

0.8 0.5

A

M

J

J

A

TKE

-40

OLD A

S

M

J

-60

S

Year: 2001

Year: 2001

-10

1.8 1.5

OBS.

1.3

TKE

1.0 OLD

0.8 A

M

J

J

A

-20

TKE

-30 -40

OLD

-50

S

-60

A

M

J

J

A

S

Month 0

2.0

Year: 2002

Year: 2002 -10

Bias (M-O)/O (%)

1.8 OBS.

1.5 1.3 1.0

TKE

-20 -30

TKE

-40

OLD

-50

0.8 0.5

A

Month

Month Concentration of NO2(µg(N)m-3)

J

0

2.0

Bias (M-O)/O (%)

Concentration of NO2(µg(N)m-3)

-30

-50

Month

0.5

-20

A

M

J

J

A OLD

S

Month

-60

A

M

J

J

A

S

Month 3

Fig. 1. The OLD versus TKE scheme. Comparison of (a) the modelled and observed NO2 in air (mg (N) m ) concentrations and (b) their biases in the period AprileSeptember for the years 1999, 2001 and 2002. M and O denotes modelled and observed value, respectively.

Future studies with this non-local scheme may focus on the interaction of the scheme with the hydrological cycle, for example, deep and shallow convection, and the land-surface parameterization.

Acknowledgement The research work described in this paper has been funded by the Serbian Ministry of Science and Environmental Protection under the project ‘‘Modelling and numerical simulations of complex physical systems’’, No. ON141035 for 2006e2010, and also by the National Science Foundation, USA. The work on this paper was partly realised by the first author during his visit to the Norwegian Meteorological Institute in Oslo.

References Alapaty, K., 2003. Development of two CBL schemes using the turbulence velocity scale. In: Fourth WRF Users’ Workshop, June 25e27, Boulder, Colorado. Alapaty, K., Alapaty, M., 2001. Development of a Diagnostic TKE Schemes for Applications in Regional and Climate Models Using MM5. Research Note. MCNC-North Carolina Supercomputing Center, Research Triangle Park, NC, p. 5. Blackadar, A.K., 1979. Modeling pollutant transfer during daytime convection. In: Fourth Symposium on Atmospheric Turbulence Diffusion and Air Quality. American Meteorological Society, Reno, NV, pp. 443e447. Berge, E., Jacobsen, H.A., 1998. A regional scale multi-layer model for the calculation of long-term transport and deposition of air-pollution in Europe. Tellus 50, 205e223. Bjorge, D., Skalin, R., 1995. PARLAM e The Parallel HIRLAM Version of DNMI. Research Report No. 27. ISSN: 0332-9879. Norwegian Meteorological Institute, Oslo, Norway. Businger, J.A., Izumi, Y., Bradley, E.F., 1971. Flux profile relationships in the atmospheric surface layer. Journal of Atmospheric Sciences 28, 181e189.

D.T. Mihailovic, K. Alapaty / Environmental Modelling & Software 22 (2007) 1685e1689 Fagerli, H., Eliassen, A., 2002. Modified parameterization of the vertical diffusion. In: Transboundary Acidification, Eutrophication and Ground Level ozone in Europe. EMEP Summary Status Report, Research Report No. 141, Norwegian Meteorological Institute, Oslo, Norway, p. 74. Hass, H., Jacobs, H.J., Memmesheimer, M., Ebel, A., Chang, J.S., 1991. Simulation a wet deposition case in Europe using European Acid Deposition Model (EURAD). In: Air Pollution Modeling and Its Applications VIII. Plenum Press, New York, pp. 205e213. Holtslag, A.A.M., Boville, B.A., 1993. Local versus nonlocal boundary layer diffusion in a global climate model. Journal of Climate 6, 1825e1842. Holtslag, A.A.M., de Bruin, E.I.F., Pan, H.-L., 1990. A high resolution air mass transformation model for short-range weather forecasting. Monthly Weather Review 118, 1561e1575. Hong, S.Y., Pan, H.L., 1996. Nonlocal boundary layer vertical diffusion in a medium-range forecast model. Monthly Weather Review 124, 2322e2339. Lenschow, D.H., Li, X.S., Zhu, C.J., 1988. Stably stratified boundary layer over the Great Plains. Part I: mean and turbulent structure. BoundaryLayer Meteorology 42, 95e121. Mihailovic, D.T., Rao, S.T., Alapaty, K., Ku, J.Y., Arsenic, I., Lalic, B., 2005. A study of the effects of subgrid-scale representation of land use on the boundary layer evolution using 1-D model. Environmental Modelling and Software 20, 705e714.

1689

Moeng, C.-H., Sullivan, P.P., 1994. A comparison of shear and buoyancy driven planetary-boundary-layer flows. Journal of Atmospheric Science 51, 999e1022. O’Brien, J.J., 1970. A note on the vertical structure of the eddy exchange coefficient in the planetary boundary layer. Journal of Atmospheric Science 27, 1213e1215. Pielke Sr., R.A., 2002. Mesoscale Meteorological Modeling, second ed. Academic Press, 676 pp. Simpson, D., Fagerli, H., Jonson, J.E., Tsyro, S., Wind, P., Tuovinen, J.-P., 2003. Transboundary Acidification, Eutrophication and Ground Level Ozone in Europe. Part I: Unified EMEP Model Description. EMEP Status Report 2003. The Norwegian Meteorological Institute, Norway, p. 74. Troen, I., Mahrt, L., 1986. A simple model of the atmospheric boundary layer; sensitivity to surface evaporation. Boundary-Layer Meteorology 37, 129e148. Zhang, D., Anthes, R.C., 1982. A high-resolution model of the planetary boundary-layer-sensitivity tests and comparisons with SESAME-79 data. Journal of Applied Meteorology 21, 1594e1609. Zhang, C., Randall, D.A., Moeng, C.-H., Branson, M., Moyer, M., Wang, Q., 1996. A surface parameterization based on vertically averaged turbulence kinetic energy. Monthly Weather Review 124, 2521e2536. Zhang, K., Mao, H., Civerolo, K., Berman, S., Ku, J.-Y., Rao, S.T., Doddridge, B., Philbrick, C.R., Clark, R., 2001. Numerical investigation of boundary layer evolution and nocturnal low-level jets: local versus non-local PBL schemes. Environmental Fluid Mechanics 1, 171e208.