boundary layer height determination in the EMEP model

Abstract. This paper introduces two changes of the turbulence parameterization for the EMEP (European Monitoring and Evaluation Programme) Eulerian air pollution model: the replacement of the Blackadar in stable and O'Brien in unstable turbulence formulations with an analytical vertical diffusion profile (K(z)) called Grisogono, and a different mixing height determination, based on a bulk Richardson number formulation (RiB). The operational or standard (STD) and proposed new parameterization for eddy diffusivity have been validated in all stability conditions against the observed daily surface nitrogen dioxide (NO2), sulphur dioxide (SO2) and sulphate (SO42−) concentrations at different EMEP stations during the year 2001. A moderate improvement in the correlation coefficient and bias for NO2 and SO2 and a slight improvement for sulphate is found for the most of the analyzed stations with the Grisogono K(z) scheme, which is recommended for further application due to its scientific and technical advantages. The newly extended approach for the mechanical eddy diffusivity is applied to the Large Eddy Simulation data focusing at the bulk properties of the neutral and stable atmospheric boundary layer. A summary and extension of the previous work on the empirical coefficients in neutral and stable conditions is provided with the recommendations to the further model development. Special emphasis is given to the representation of the ABL in order to capture the vertical transport and dispersion of the atmospheric air pollution. Two different schemes for the ABL height determination are evaluated against the radiosounding data in January and July 2001, and against the data from the Cabauw tower, the Netherlands, for the same year. The validation of the ABL parameterizations has shown that the EMEP model is able to reproduce spatial and temporal mixing height variability. Improvements are identified especially in stable conditions with the new ABL height scheme based on the RiB number.


Introduction
Air quality models are nowadays recognized as an important tool for air quality assessment.Although measurements are the basis of air quality assessment, there are several advantages provided by models: high spatial and temporal resolution of simulated data, forecasting of air quality as a result of changes in emissions or/and meteorological conditions and a better understanding of the physical processes that drive the transport of pollutants in the atmosphere.For nearly 30 years, the European Monitoring and Evaluation Programme (EMEP) under the Convention on Long-Range Transboundary Air Pollution (LRTAP), has been responsible for the development of air quality modelling systems to support the design of the environmental control strategies in Europe.The Unified EMEP model (Simpson et al., 2003) was developed and used to simulate transboundary transport of air pollution on the European scale.Recently, special applications of the model have been developed at higher resolutions, and coupled with different meteorological drivers: EMEP4UK (EMEP for the United Kingdom; e.g.Vieno et al., 2009;Vieno et al., 2009) and EMEP4HR (EMEP for Croatia; e.g.Jeričević et al., 2007;Kraljević et al., 2008).Development of the EMEP model includes detailed meteorological effects that become progressively more important on the finer spatial scale, such as turbulence and convection Published by Copernicus Publications on behalf of the European Geosciences Union.generated by a complex terrain.Turbulence parameterizations, particularly schemes for calculation of vertical diffusion coefficients, K(z), where z is the height, needs to be tested as a first step of the EMEP model development on a finer horizontal scale.
Previous studies have already shown that the parameterizations of K(z) have significant impacts on the simulated chemical concentrations (e.g.Nowacki et al., 1996;Biswas and Rao, 2000;Olivie et al., 2004).Different parameterizations for K(z), depending on the stability in the atmospheric boundary layer (ABL), have been proposed (e.g.O 'Brien, 1970;Deardorff, 1972;Louis, 1979;Holtslag and Moeng, 1991;Holtslag and Boville, 1993;Grisogono, 1995).O'Brien (1970) suggested a simple parameterization K(z) scheme used in many air quality models ranging from simple 1-D models (e.g. Lee and Larsen, 1997) towards application as in complex chemical models e.g.Comprehensive Air Quality Model with Extensions (CAMx, http://www.camx.com/; ENVIRON, 1998;Zhang et al., 2004), as well as in the EMEP model (Fagerli and Eliassen, 2002).In CAMx there are a few K(z) parameterization schemes, with the O'Brien scheme as one of the options.Presently, in the EMEP model the O'Brien scheme is used for the convective boundary layer (CBL), while in the stable boundary layer (SBL) conditions Blackadar (1979) scheme that is based on Monin-Obukhov (M-O; Monin and Obukhov, 1954) similarity theory is applied.There are many studies which show that the surface-layer formulations based on the M-O theory are often not applicable in statically stable conditions (e.g.Mahrt, 1999;Pahlow et al., 2001;Poulos and Burns, 2003;Mauritsen et al., 2007;Grisogono et al., 2007).A new proposed non-local scheme, called Grisogono, is implemented in the model and it is not based on the M-O similarity theory.The Grisogono scheme, applied in SBL and CBL, uses a linear-exponentially decaying profile, generalizing the O'Brien third-order polynomial K(z) (Grisogono andOerlemans, 2001a, b, 2002).The Grisogono method uses empirical coefficients determined on the LES data in stable and neutral conditions and presently the same coefficients are employed in the model for convective conditions.Previous work from Jeričević and Večenaj (2009;JV09) is now extended with a new integral empirical coefficients which are more convenient for application in numerical and air quality models.
Special emphasis is given to the ability of the ABL height scheme to capture the vertical transport and dispersion of atmospheric air pollution.A significant influence of the ABL height (H ) on the surface nitrogen oxide (NO x ) and the particulate matter (PM) concentrations has often been found in urban and suburban areas (e.g.Schäfer et al., 2006), while Athanassiadis et al. (2002) show that an accurate H determination is needed to properly simulate pollutant levels with grid-based photochemical models.Furthermore, H is explicitly included in the both EMEP K(z) parameterizations.Therefore it is important to evaluate the EMEP model ability to simulate the spatial and temporal variability of H .The operational (e.g.Jakobsen et al., 1995;Seibert et al., 2000) and a new ABL height scheme based on the bulk Richardson number (Ri B ) are evaluated.The Ri B method is a standard and widely used approach to derive H from the numerical weather prediction (NWP) models, as well as from the radiosounding data (e.g.Mahrt, 1981;Troen and Mahrt, 1986;Sørensen et al., 1996;Fay et al., 1997;Seibert et al., 2000;Zilitinkevich and Calanca, 2000;Zilitinkevich and Baklanov, 2002;Gryning and Batchvarova, 2002;Jeričević and Grisogono, 2006).
The operational version of the EMEP model, and the version with new parameterization schemes (i.e.K(z) and ABL height schemes) are verified by comparing one full year of the modelled data against the corresponding set of measurements from different EMEP stations in Europe.Based on this validation, discrepancies (both in the measurements and in the model) are identified.Pronounced differences between the performances of the two model versions and impacts on the simulated concentrations are investigated and recommendations for future work are provided.This paper gives the basis for further development and improvement of the EMEP model by e.g.improving the parameterizations of the vertical diffusion and the boundary layer representation.This study has been conducted within the EMEP4HR project whose main purpose is to develop and test an operative framework for the environmental control of air pollution problems in a broader region of Croatia.Previous efforts addressing the same issue are described in Klaić (1990Klaić ( , 1995Klaić ( , 2003)), and Klaić and Beširević (1998).To summarize, we try to combine several recent findings about the nature, theory and modeling of the ABL in an operational atmospheric chemistry model.

The EMEP model description
The Unified EMEP model (http://www.emep.int/)was developed at the Norwegian Meteorological Institute under the EMEP programme.The model is a development of the earlier EMEP models (Berge and Jakobsen, 1998;Jonson et al., 1998), and is fully documented in Simpson et al. (2003) and Fagerli et al. (2004).The model has been extensively validated against measurements (Fagerli et al., 2003(Fagerli et al., , 2007;;Simpson et al., 2006aSimpson et al., , b, 2007;;Jonson et al., 2006;Tsyro et al., 2007;Fagerli and Aas, 2008).It simulates the atmospheric transport and deposition of acidifying and eutrophying compounds, as well as photo-oxidants and particulate matter over Europe.The model domain covers Europe and a part of the Atlantic Ocean with the grid size 50 km×50 km while in the vertical there are 20 terrain following layers reaching up to 100 hPa.The Unified EMEP model uses the 3-hourly meteorological data from the PARallel Limited

Measurements
Different data sets have been used here to evaluate the EMEP model performance: (i) observed daily surface concentrations of NO 2 , SO 2 and SO 2− 4 at different EMEP stations in Europe during the year 2001 (Fig. 1), (ii) radiosounding measurements from various European cities in January and July 2001 (Table 1) and (iii) wind and temperature profiles from the Cabauw tower, the Netherlands, also in the year 2001.
The selected pollutants are among the most important acidifying and eutrophying pollutants contributing to air pollution and atmospheric chemistry.Sulphate is a secondary pollutant, an oxidation-product of SO 2 , which contributes to acid rain formation.Since atmospheric lifetimes of SO 2 and NO 2 are 1 to 3 days and their oxidation product's lifetime is generally even longer (Seinfeld and Pandis, 1998), they are subjected to the atmospheric transport and mixing processes, and therefore suitable for validation of vertical diffusion scheme efficiency.Furthermore, NO 2 , SO 2 and SO 2− are monitored at the majority of EMEP stations with a good spatial and time resolution.

Measurements from the EMEP stations
This study has used the measurements at the EMEP stations (http://www.emep.int/)for the model evaluation.These measurements are well documented, quality controlled and they mostly represent background conditions over a larger area.In order to obtain data that are characteristic for long-range transport, it is important that a station is representative of the EMEP 50×50 km 2 grid square averages.It should be emphasised that the recommendation for the EMEP sites not to be influenced by local pollution implies that their location is chosen to ensure representativeness of the lower concentrations in the grid, not the grid average.Also, the measurements are not of equal quality at all stations, and to some extent, this fact may be explained by different measurement methods (e.g.Fagerli et al., 2003).
The analyzed stations within the EMEP domain are shown in Fig. 1.Most of the stations are below 300 m (blue dots).Nevertheless, many stations in the Central European area are located between 600 m and 1000 m, while in the Alps area stations are often above 1000 m.Jungfraujoch (CH01) in Switzerland is above 3000 m and Chopok (SK02) in Slovakia is above 2000 m and they are not used for the evaluation of turbulence parameterization schemes.Mountain stations are not very well represented in models with coarse horizontal resolution, having too low altitude in the model and consequently, surface concentrations are too high compared to measurements.The orography misrepresentation is a known modelling problem (e.g.Žagar and Rakovec, 1999;Ivatek Šahdan and Tudor, 2004).
A list of all EMEP stations with more details on the measuring programme and available data can be found at: http:// tarantula.nilu.no/projects/ccc/network/index.html.The number of used stations varied from compound to compound i.e. the measured daily SO 2 was available at 68 stations, NO 2 at 43 stations and SO 2− 4 at 58 stations.

Measurements from the radiosounding stations
Radiosoundings are often used in order to operationally determine and verify H values (e.g.Seibert et al., 2000).Nevertheless, these measurements are usually only taken twice a day at 00:00 and 12:00 UTC and consequently the soundings can only be used as an overall reference.The data possess reasonably good spatial distribution over Europe and they are commonly available and quality controlled.In this paper, the evaluation was performed using the data obtained from 24 different measuring stations in Europe (Table 1) during January and July in 2001.

Cabauw measurements
The Cabauw tower is located in the western part of the Netherlands (51 • 58 N, 4 • 56 E) with flat surroundings e.g.van Ulden and Wieringa (1996).Temperature and wind averages are computed over 10 -min intervals.Wind speed and wind direction are measured at six levels: 10, 20, 40, 80, 140 and 200 m while temperature is measured at one additional level, at 1.5 m.Pressure is measured at 1.5 m height only.A hydrostatic balance is assumed in order to derive potential temperature needed for the Ri B .Pressure on upper levels is integrated from the surface pressure at 1.5 m using the trapezoidal rule.The Cabauw observations have been used in other studies to validate the land surface parameterization schemes (e.g.Beljaars and Bosveld, 1997;Chen et al., 1997;Ek and Holtslag, 2005).
The measurements from the Cabauw tower have a high resolution in time and their vertical distribution is dense enough to reconstruct physical processes in the surface layer (occasionally even higher) thus providing the possibility to investigate and analyze the ABL structure near the surface into greater detail than with "standard" measurements.

The LES data
Data from the DATABASE64 (Esau and Zilitinkevich, 2006) is used in order to evaluate the performance of different K(z) schemes, in stable and neutral atmospheric conditions and to determine empirical coefficients applied in the Grisogono approach.The LES used a dynamic sub-grid scale closure model which parameterizes turbulent kinetic energy (TKE) dissipation with Smagorinsky closure and a resolution of 64 3 gridpoints.Esau and Zilitinkevich (2006) show on a few intercomparison and convergence studies, that relatively small 64 3 mesh is sufficient to keep simulation errors at the level less than 5% of the total turbulent kinetic energy.Furthermore, the authors have provided comparisons for some turbulence statistics resolved on 64 3 and 128 3 meshes finding that the differences in the vertical transport characteristics remains fairly small (e.g. for weakly or moderately stable ABLs).They found this conclusion consistent with the Beare et al. (2006).
The DATABASE64 consists of a wide range of neutral and stably stratified cases.In all cases the initial temperature profile (neutral or with constant stratification), the constant background geostrophic wind, the surface roughness length and surface heat flux were defined.It should be pointed out that Basu et al. (2008) showed, based on an analytical approach, that application of the surface heat flux should be avoided as a lower boundary condition in LES models due to the existence of "dual" nature of sensible heat flux in stable conditions.Different stability classes are defined in the DATABASE64: CNT -conventionally neutral cases; TSTtruly stable cases i.e. nocturnal stable; CST -conventionally stable cases i.e. long-lived stable cases.

The determination of integral empirical coefficients
Since we deal here with 3-D realistic flows, transport and dispersion using the EMEP model, and thus departing from certain idealizations in the LES results (Esau and Zilitinkevich, 2006;Basu et al., 2008;JV09), we conveniently extend and generalize our estimates of eddy diffusivity for momentum (K m ) in stable and neutral conditions compared to that in JV09.This generalization for the turbulent momentum Table 2.The new integral empirical coefficients for maximum of K(z) and its height, C(K) and C(z max ) respectively, both for momentum and heat turbulent transport, with the corresponding standard deviations σ , used for determination of K(z) profiles with the Grisogono approach.
transfer goes from the following equation (applicable to nearsurface turbulent processes only) which was used in JV09: where k is the von Karman constant, k≈0.41, u * is a friction velocity (m s −1 ).Now Eq. ( 1) is commonly extended upward with: where U is the mean wind speed and the u w is the vertical sub-grid scale flux of the momentum.The eddy diffusion for heat (K h ), remains in the flux form as in JV09: where w θ is the related vertical flux of heat.Note that Eqs.
(2) and ( 2) must be valid throughout the ABL, provided the K-theory holds true.
The old empirical coefficients for K m (Table 3 in JV09) are now extended with a new integral empirical coefficients for the maximum of K(z) and its height, (C(K) and C(z max ), respectively), which are more convenient for application in numerical and air quality models (Table 2) simply because these coefficients are now based on more generalized ABL equation (Eq. 2 instead of Eq. 1), and more properly chosen LES data subsets.These should, we believe, better pertain to the overall needs, concept and resolution of the EMEP model.Standard deviations of the new coefficients are also provided.The new C(z max ) is simply an inverse of the previous C(h) value in JV09.
The earlier results (JV09) showed greater efficiency of the momentum transport Eq. ( 1), relative to the heat transport Eq. ( 3).The integral empirical coefficients for momentum differ significantly from that in JV09 mainly due to application of two different equations for K m (Eqs. 1 and 2).The main result here is that with the new integral empirical coefficients (Table 2), a better balance between the momentum and sensible heat is achieved.For heat the same equation, i.e.Eq. ( 3), is used in JV09 and in this work.However, the coefficients C(K) for heat (Table 2) differs slightly from the old ones (Table 3 in JV09), because a few LES runs with a higher uncertainty are omitted in the new calculations.
In this work statistical evaluation of the LES runs prior to the calculation of the empirical coefficients is provided while in JV09 the only criterion on set of the LES data is the size of the domain.The variability intervals of the new empirical coefficients are shown in Fig. 2. Obviously, the variability of the empirical coefficients for the mechanical turbulent transport is within 50% of the average value while for the heat transport, estimated variability is somewhat lower, between 38% and 40%.

The operational K(z) schemes in the EMEP model
In the EMEP model above the ABL and inside the SBL K(z) is calculated with the local scheme proposed by Blackadar (1979): where l is the turbulent mixing length (m), V H is horizontal wind speed, z is the model layer thickness, | V H / z| is the absolute value of wind shear in the vertical.In the EMEP model l is parameterized crudely according to: where z is the height above the ground and z m =200 m.As a side note with l from Eq. ( 5), it is notoriously difficult to parameterize stratified flows, and l is invariably treated poorly in modelling the SBL (e.g.Grisogono and Belušić, 2008;Grisogono, 2010).The R i is the gradient Richardson number defined as: where n is the model level, θ is a potential temperature (K), θ is the potential temperature difference in the model layer, and Ri C is the critical Richardson number calculated from McNider and Pielke (1981): where A=0.115, B=0.175 and z 0 =0.01 m.The final Ri C value is: Ri C = MAX 0.25,0.115(z) 0.175 .Obviously with z→0, Ri C →0.25.
In the unstable ABL, K(z) is calculated with the O'Brien scheme: where K H is a K(z) value at the top of the ABL, i.e.K(z=H ) and K H S is a K(z) value at the top of the surface-layer (H S ).It is assumed that ∂K(z)/∂z=0 at z=H .From the M-O similarity theory for the surface layer (e.g.Stull, 1988): where L is the Monin-Obukhov length (m) and is a universal function.The friction velocity is given by: where τ is the near-surface turbulent momentum flux (N m −2 ) and ρ (kg m −3 ) is air density (derived from surface pressure and temperature).The L is given by: where θ S is a surface potential temperature, Q h is the sensible heat flux (W m −2 ) taken from the NWP PARLAM-PS model, g is acceleration due to gravity (9.8 m s −1 ) and C p is a specific heat capacity of dry air at constant pressure (1005 J kg −1 K −1 ).Universal functions used in the EMEP are those recommended by Garratt (1992) for statically unstable cases: and for statically stable cases: To avoid the non physically small exchange coefficients within the ABL, K(z) value is evaluated at the top of the lowest model layer (z≈90 m) with Eq. ( 9) both in stable and unstable atmospheric conditions.

Definition of the Grisogono K(z) scheme
A gradually-varying function, i.e.K(z), was introduced to generalize the classical analytical solution for the Ekman layer flow using the WKB method (after Wentzel, Kramers and Brillouin who popularized this method in theoretical physics).Since there is no explicit relation between the ABL profiles and K(z), a solution which generalizes the third order O'Brien polynomial was defined between a constant K value and a numerically derived solution of the Ekman profile (Grisogono, 1995).Furthermore, the Prandtl model for katabatic flows is solved for gradually varying K(z) expressed in an exponentially decaying form (Grisogono andOerlemans, 2001a, b and2002).This newly proposed scheme where the O'Brien thirdorder polynomial K(z) is generalized into a linearexponentially decaying function (e.g.Grisogono and Oerlemans, 2002):  where z max is the height of K(z) maximum value (K max ).The Grisogono profile combines a linear term, which dominates near the surface, with an exponential decay, so that the maximum of K(z) is reached at about 0.3 H , similar to the O'Brien's formula.One can notice that one of the advantages of Eq. ( 12) over the O'Brien's Eq. ( 8) is that it needs only two input parameters, K max and z max .Parameters z max and K max are evaluated from the following equations: The equation for the height z max of K max (Eq.8 in JV09) is slightly modified in order to get a linear dependence between the z max and C(z max ).The new expression is: where C(K)=0.05 and C(z max )=0.21 are the new integral empirical constants for K h estimated from the LES data (Table 2).By inserting Eqs. ( 13) and ( 14) into ( 12), a new simplified form is derived: where C≈0.39 is a new compiled constant.This new simplified form explicitly includes u * and H , utilized from the meteorological driver and its accuracy is constrained with the NWP model performance.

The ABL height
The ABL height is an important parameter, which limits the modelled vertical extent of continuous turbulent mixing in the atmosphere starting from the surface.The operational method for the calculation of H in the EMEP model determines H from the NWP PARLAM-PS output (Jakobsen et al., 1995;Simpson et al., 2003).In statically stable conditions H is calculated as the height where K(z)<1 m 2 s −1 , with K(z) profiles calculated with the local Blackadar method, Eq. ( 4), and vertically linearly smoothed over a few adjacent layers.As mentioned in the introduction, the Ri B method is a standard and widely used approach to derive H from the numerical weather prediction (NWP) models, as well as from the radiosounding data.The proposed and commonly used Ri B method is based on the assumption that continuous turbulence vanishes beyond Ri BC , some previously defined critical value of Ri B .The height at which Ri B reaches Ri BC is considered as H .It is defined as: Here θ 1 is a potential temperature at the lowest model level, z 1 , and θ (z) is an average potential temperature between heights z and z 1 .Now, H is the height where Ri BC =0.25 is reached.However, the supposed existence of Ri BC has recently been criticised (Zilitinkevich and Baklanov, 2002;Jeričević and Grisogono, 2006;Mauritsen et al., 2007;Zilitinkevich et al., 2008;Grisogono and Belušić, 2008) and the development of improved K(z) schemes based on higher order closures is a subject of current and future research.The main advantages of the new method are that Ri B includes the two major turbulence generators in the atmosphere, thermal and mechanical sources of turbulence, and it is applicable in statically stable and unstable atmospheric conditions.Equation ( 16) describes H as an integral atmospheric property that relates surface processes to the upper-level processes in the ABL and thus comprises non-local effects.The main weakness of the operational ABL height method in stable conditions is the dependence on the K(z) profiles calculated with the Blackadar approach Eq. ( 4).The operational method in statically stable conditions is based on the Ri number and it also includes both sources of turbulence; however, it can be oversensitive to the local turbulence and may underestimate the ABL height.In statically unstable conditions, the accuracy of the operational method depends on surface parameters obtained from the NWP model e.g.Q h , and vertical distribution of Q h via dry adiabatic adjustment, while effects of the mean wind shear are not included.

Statistical methods
The correlation coefficient (r) and BIAS= M−O O ×100%, are calculated between the observed (O), daily surface NO 2 , SO 2 and SO 2− 4 concentrations (c(NO 2 ), c(SO 2 ), c(SO 2− 4 )), and the corresponding modelled values (M).Furthermore differences (D) between the old and new r and BIAS values are calculated in order to find potential improvements in the EMEP model performance with the change of K(z) scheme in all stability conditions.Differences are defined as: where parameter X can be r or the absolute value of BIAS, (ABS(BIAS)).For X=r, D(r)>0 means that the model performs better with the Grisogono K(z) scheme, while for X=BIAS, D(BIAS)>0 denotes that the STD scheme agrees better with the observations.Similarly D≈0±0.001denotes an equally good performance of the both schemes.The modelled absolute values and BIAS are very sensitive to the balance between the different processes in the model.Therefore, a smaller BIAS between the model and measurements does not necessarily mean that the new scheme is better than the standard one; it only means that the average concentrations determined with the new scheme are closer to the average of the observed concentrations.However, the BIAS can give an insight into the general effect of the new scheme on the modelled values.For instance, if the Grisogono parameterization is less diffusive in statically stable conditions this should lead to higher average concentrations in these cases.The temporal correlation coefficient, however, is a better measure for whether the new scheme provides a better physical description of overall simulated concentrations.Therefore, we focus on the changes in the correlation coefficient between the model results and observations.

Significance tests
A standard significance Fishers z-test (Fischer, 1915) is conducted in order to find whether the changes in r, which resulted due to variations of the K(z) and ABL height schemes in the EMEP model, reflect the change of stochastic relation between the two data sets.
The hypothesis H 0 : r 1 =r 2 , and H 1 :r 1 =r 2 , have been tested, where r 1 and r 2 are the correlation coefficients determined between the observations and modelled data calculated with two different K(z) schemes, the STD and Grisogono and with the two different ABL height schemes.For the 95% confidence interval hypothesis ≤2 is satisfied.Variables z 1 , z 2 and σ z 1 −z 2 are determined from: where n 1 and n 2 are the sizes of analyzed data sets.However, the appropriateness of this procedure is questioned since initial assumptions for its application are not completely satisfied, i.e. the mutual independence of the observation and modelled data, and the distribution of the quantity following a normal distribution.The z-test has been used in practice, nevertheless it is found to be quite insensitive to establishing whether two correlations have different strengths.In this test, as in many other standard statistical tests, an assumption of mutual independence is made.However, daily concentrations are not completely independent since they are time-correlated with the persistence of Fig. 3.The average vertical diffusion profiles of momentum K Gm (solid black line) and heat K Gh (solid red line) calculated with the Grisogono method and the average K(z) profiles calculated with the O'Brien method K OB (dashed gray line) against the average vertical profiles of eddy diffusivity for momentum K m (black dots) and heat K h (red dots) estimated from the LES data.The LES runs are averaged according to stability (CNT -conventionally neutral cases; TST -truly stable cases i.e. nocturnal stable; CST -conventionally stable cases i.e. long-lived stable cases) and the ABL height, H , i.e.: H <300 m, 300<H <600 m, and H >600 m. meteorological events (Fox, 1980;Chang and Hanna, 2003).Time correlation in data sets may affect significance tests in many different ways making the estimation of degrees of freedom needed for the level of significance determination impossible.Willmott (1982) argued that it is inappropriate to report r as statistically significant, among other reasons because the magnitude of r and its associated significance level are not necessarily related to the accuracy of the simulated concentrations, and rarely conform to the assumptions that are prerequisite to the appropriate application of inferential statistics, as it was also stated here.

The K(z) profiles from the LES data
New K(z) profiles, in stable and neutral conditions, are calculated using the empirical coefficients for K m and K h from Table 2.The average K m and K h profiles calculated with the Grisogono and O'Brien methods are plotted against the corresponding profiles estimated from the LES data (Fig. 3).The O'Brien and Grisogono schemes are compared here because they are similar non-local methods, although the O'Brien is not applied in the stable conditions in the EMEP model.The LES runs are averaged according to stability (CNT, TST and CST stability classes) and ABL height, i. -medium, with 300 m<H <600 m, and -high, with H > 600 m.
The stability and H classes are marked above every profile in Fig. 3.The Grisogono is in a good agreement with the LES in the CNT conditions, especially for H <600 m, while for H >600 m the Grisogono method underestimates the turbulent transport for the heat and momentum.The O'Brien method overestimates momentum and heat transport compared to the LES data for all H classes during the CNT neutral conditions (Fig. 3).The agreement between the Grisogono and O'Brien for momentum transfer is good, especially for H <600 m during the nocturnal stable conditions (TST).However, both schemes overestimate K m .The O'Brien scheme efficiently represents heat transport during the TST conditions for the H >600 m, while the Grisogono agrees better with the K m .The strongest stability occurs in the long lived stable class i.e.CST where H is mainly <300 m and the average eddy diffusion does not exceed 1.5 m 2 s −1 .The O'Brien and Grisogono perform similarly slightly overestimating the K m and K h in the CST.
Since the O'Brien scheme is used in the EMEP only in convective conditions the focus here is mainly given to the intercomparison of the Grisogono and Blackadar approach which is applied in the EMEP in stable conditions.Therefore, the Blackadar method, Eq. ( 4), is applied to the LES data and K(z) profiles are determined (K Black ) and compared to the LES, O'Brien and Grisogono profiles (Fig. 4).Firstly and obviously, the Blackadar method severely overestimates K max for most of the cases.However, note a relatively good agreement between the K Black and K h >30 m 2 s −1 in neutral and moderately stable LES cases with H >600 m, though there are fewer cases there.It is known that local schemes, such as Blackadar, describe well phenomena like residual layer, low level jet or clear air turbulence above the ABL (Stull, 1988).The simulated nocturnal SBL develops in a near neutral background atmosphere, with heat loss at the surface, and occurs during night time over land with a nearneutral residual layer.The area of the intensified local mixing can be a residual of the convective mixing or a low level jet resulted from wave breaking, or other intensive forces.

The evaluation of the K(z) schemes performance in
the EMEP model

The EMEP network data
In order to quantify changes in the model performance in all stability conditions acquired with the new K(z) scheme, differences between the correlation coefficients, D(r) in Eq. ( 19), obtained by two different parameterization schemes are calculated for NO 2 , SO 2 and SO 2− 4 .The spatially interpolated annual correlation coefficients (a) r(NO 2 ), (b) r(SO 2 ) and (c) r(SO 2− 4 ) for the operational EMEP model, and the spatially interpolated differences in the annual rvalues, D(r), acquired with the new K(z) scheme (d) D(r(NO 2 )), (e) D(r(SO 2 )) and (f) D(r(SO 2− 4 )) are shown in Fig. 5.The upper panels in Fig. 5 show the operational model performance according to r values.The model is in a good agreement with the measurements finding correlation coefficient 0.5≤r(NO 2 )≤0.75 at 56% stations, while at 44% stations r(N O 2 )≤0.5 (Fig. 5a).For SO 2 it is 0.5≤r(SO 2 )≤0.77 at 43% stations, and r(SO 2 )≤0.5 on 57% (Fig. 5b).For SO 2− 4 it is 0.5≤r(SO 2− 4 )≤0.87 at 86% stations while only at 14% of the analyzed stations r(SO 2− 4 )≤0.5 (Fig. 5c).It should be pointed out that r(SO 2− 4 ) is the highest among all analyzed species with r(SO 2− 4 )>0.7 at 31% stations.The operational EMEP model performance has been regularly evaluated by comparison with observations of air and precipitation data compiled in the EMEP network.Results of the model evaluations have been published in the official reports (http://www.emep.int/publications.html).The analyzed year was not exceptional regarding meteorological conditions and the EMEP model performance is in agreement with the previous evaluation results (Fagerli et al., 2003).
The lower panels in Fig. 5 show improvements (blue colour) and deteriorations (red colour) in r values as a consequence of different K(z) scheme employment in the EMEP model.Although generally an improvement is detected, there are still some areas where the STD method has a better performance.Better results with the STD scheme are found for NO 2 in Scandinavian area and Italy; likewise for SO 2 , only for the stations in the northern part of Great Britain.For the sulphate similar or lower results with the Grisogono scheme are obtained in Scandinavia, Great Britain and Hungary.However, the spatial interpolation analysis should be carefully interpreted since the results may be influenced by a few stations in the areas with low resolution in the measurements (here Central and Eastern Europe).

Yearly analyzes at selected stations
As previously explained in Sect.2.3.1, the EMEP recommends that measuring stations are located away from a large local emission sources.If the station is affected by local sources, irregular variability is observed in concentrations, which is not modelled with the EMEP model, and underestimation as well as overestimation of the measurements may occur.
Based on the operational EMEP model evaluation in the year 2001 discrepancies between the model and measurements are identified.Discrepancies with a factor of two or more between the model and measurements are found at different stations which can be categorized as: (i) stations where peak events or episodes occurred in the measurements influenced by local emission sources, and stations in the vicinity of large emission sources (e.g.shipping area in the North Sea) and (ii) mountain stations.Since shipping emission paths are not sufficiently resolved due to the coarse horizontal resolution in the model, higher concentrations are horizontally diffused over larger areas (including analyzed stations, where obviously these high concentrations were not observed).Generally, stations in the North Sea shipping area are probably overestimated with the EMEP model due to the coarse model horizontal resolution, but it might be due to other reasons e.g.emissions, meteorology, chemistry, etc. Stations with the highest discrepancies were excluded from the annual r and BIAS estimation.
The NO 2 time series are analyzed all stations except those with higher discrepancies in order to investigate seasonal variability of K(z) with the two different schemes applied.The annual course of (a) r values, (b) BIAS values, (c) RMSE and (d) average monthly concentrations of NO 2 calculated between the measurements and modelled c(NO 2 ) values with two K(z) schemes, the Grisogono (blue line) and STD (red line) are shown in Fig. 6.In Fig. 6a systematically higher r values with the new K(z) scheme are shown in both: statically stable conditions, more characteristic during the colder part of the year, and statically unstable conditions, during the warmer part of the year.According to BIAS (Fig. 6b), in the warmer part of the year the model underestimates c(NO 2 ) with both K(z) schemes.Furthermore, RMSE in Fig. 6c is also the lowest during the summer time.The measured and modelled mean monthly NO 2 values in Fig. 6d show a decrease of c(NO 2 ) during the warmer part of the year.This drop in c(NO 2 ) is caused by the increased photolysis of NO 2 and more vigorous vertical mixing during the warmer period.A seasonal variation of NO 2 emissions also plays a significant role in the annual c(NO 2 ) course.Note a higher c(NO 2 ) values with the new K(z) scheme during the warmer part of the year, which shows that the new K(z) scheme is less diffusive in convective conditions than the STD scheme.In Fig. 6d note that average monthly values with the both schemes are similar during the colder part of the year, while the second peak in November is not captured with the model.Nevertheless, r is higher with the new scheme in stable winter conditions as well.
The annual r and BIAS values between the measured and modelled daily surface c(NO 2 ), c(SO 2 ) and c(SO 2− 4 ) concentrations are also calculated (not shown).With the Grisogono scheme r(NO 2 )=0.65, r(SO 2 )=0.57, while r(NO 2 )=0.63, r(SO 2 )=0.55 are attained with the STD method.For sulphate both schemes have a similar result, r(SO 2− 4 )≈0.64.According to the BIAS values the model generally overestimates SO 2 ≈27% with the Grisogono and ≈30% with the STD method.The model underestimates SO 2− 4 and NO 2 , BIAS (SO 2− 4 )≈−19% with the original scheme and BIAS (SO 2− 4 )≈−13% with the new K(z) scheme, while BIAS (NO 2 )≈−18% with both schemes.The overestimation of SO 2 and the underestimation of sulphate indicate that other processes responsible for sulphate formation in the model Atmos.Chem. Phys., 10, 341-364, 2010 www.atmos-chem-phys.net/10/341/2010/should be investigated as well as meteorology, particularly precipitation and moisture provided by the NWP model.These evaluation results suggest that the new scheme is at least marginally better, and definitely simpler, than the STD scheme.

The ABL height representation in the EMEP model
The operational and the new ABL height scheme based on the Ri B number are compared in the EMEP model.The evaluation is performed on two data sets: (i) radiosoundings from 24 different measuring stations in Europe (Table 1) during January and July in the year 2001 and (ii) on vertical temperature and wind measurements in the year 2001 from the Cabauw tower.

Radiosounding data
The r and BIAS values are calculated between the H determined from the soundings (H sond ), and H calculated from the EMEP model (H EMEP ) for January at 00:00 UTC and July at 12:00 UTC in the 2001.The H EMEP is determined with the ABL height scheme (H std ) and also with the Ri B scheme (H new ).The values of H sond are determined with the Ri B scheme. Figure 7a shows r in January.Generally, r≈0.7 with lower values r<0.3 are found at Torshaven, Legionowo, Payerne and Izmir stations, while higher values r>0.7 are found at: Uccle, Herstmonceux, Trappes, and Stavanger.The H new shows a moderate improvement in r, while there is a considerable improvement in BIAS values, see Fig. 7b.The model underestimates H sond with the standard scheme (BIAS≈−50%), while with the new ABL height scheme the underestimation is generally significantly lower (BIAS≈−30%).Figure 7c shows the average monthly H at 00:00 UTC calculated from the soundings, H sond , with values 100 m<H sond <1000 m.The highest H sond at 00:00 UTC are found at the stations located in Southern Europe e.g.La Coruna and Lisbon.However, Torshaven, Trappes and Stavanger have higher H sond than other northern stations.On the other hand, the lowest H sond in January are found for the stations in Central Europe e.g.Payerne, Meiningen, Prague, Vienna, Wroclaw and Milan, which is expected due to long stable conditions occurring over the continent during the winter, and the corresponding H are usually low.The average H calculated from the model with the standard (H std ) scheme generally underestimates H sond (see Fig.     1) in January 2001 at 00:00 UTC.
exception is Payerne where both ABL schemes overestimated H std .Payerne is situated in the Alps and obviously the model did not manage to simulate the SBL in a complex orography.
The time series of H values in January are shown in Fig. 8 for four selected locations; two with a higher r, i.e.Herstmonceux and Stavanger (Fig. 8a and b, respectively), and two with a lower r, i.e.Torshaven and Legionowo (Fig. 8d  and c, respectively).For Herstmonceux and Stavanger the agreement between the H sond and H EMEP is good, especially with the new ABL height scheme.Note a period of low H EMEP ≈100 m (Fig. 8b, c and d   For that period at Torshaven the difference between H sond and H EMEP is ≈1000 m, while at Legionowo ≈500 m.This disagreement between H sond and H EMEP at Torshaven and Legionowo during stable conditions is the main cause for the correspondingly lower r values.
underestimations, are achieved with the standard method.
The measured average values are much higher in July at 12:00 UTC than in January at 00:00 UTC, as expected, with 1000 m<H sond <2000 m.However, the model average values are much lower for both methods varying between 400 m and 1500 m.The time series in July (Fig. 10

The ABL at the Cabauw
In this section the average hourly vertical profiles of the Ri B number ((Ri B (z j ,t))), where j=10, 20,. . ., 200 m are the measuring levels; and the corresponding H at the Cabauw tower are analyzed and described for every month in the year 2001 (Fig. 11).
As mentioned H from the Cabauw data (H tower ) is determined with the Ri B method.Vertical profiles of the Ri B number are calculated from the temperature and the wind measured at every tower level with the time interval t=10 min.In this way, the sequence of Ri B (z,t) values for the year 2001 is produced and monthly averaged to obtain the Ri B daily courses of Ri B (z j ,t) (Fig. 11).It is relatively easy to follow daily and seasonal variations of H by looking at the Ri BC =0.25 (thick blue line at the top of the blue area in Fig. 15).
The analysis of Ri B (z j ,t) provides a good insight in the processes of development and decay of the CBL and SBL at different times of the year.The occurrence of the morning and the afternoon transition layer, characterized with a sudden and rapid decay/increase of the CBL, is also shown.In January, Fig. 11a, during the nighttime, H is often less than 100 m.Daily development of H starts after 10:00 a.m.reaching the maximum H ≈200 m at 01:00 p.m. and lasting approximately 1 h after which H decreases.In February, Fig. 11b, the nighttime H is higher than in January, ranging between 100 and 200 m; the CBL starts to develop around 08:00 a.m.reaching the maximum between the noon and 02:00 p.m.In February the afternoon transition layer occurs around 03:00 p.m.Note that the transition layer has similar characteristics for the most of the analyzed months in the year 2001.In following spring and summer months from March,Fig. 11c,to August,Fig. 11h, the CBL is progressively intensifying, becoming more and more unstable.In the warmer part of the year the CBL lasts longer, which is expected since the CBL is correlated with the incoming solar radiation.Note appearance of the areas with Ri B (z j ,t)<0 numbers (from yellow to red areas in Fig. 11) in April, becoming largest in June, Fig. 11f.During the SBL conditions, in the warmer part of the year, strong near surface inversions and weak winds are measured in the surface layer.In the nigh-time SBL conditions, Ri B (z j ,t)>>Ri B C (white areas in Fig. 11) is found and the corresponding H is extremely shallow.Stable conditions prevail in September and October and SBL is 100 m-150 m thick.Dominantly stable conditions with mostly Ri B (z,t)>0 are present in November and December, Fig. 11k and l respectively.Unstable conditions occur from 10:00 a.m. to 14:00 p.m. and the average H is only 100 m.
Figure 12 shows monthly correlation coefficients calculated between the H determined from the measurements, H tower , and the modelled values determined with the operational and Ri B number method; H std (red) and H new (blue), respectively.Obviously the new ABL height scheme gives significantly better results for all months, except for February and June when both schemes performed similarly.
From Fig. 11 it is obvious that an estimated H exceeds 200 m often, especially during the warmer part of the year, which significantly limits the possibilities for the model evaluation.In Fig. S1: http://www.atmos-chem-phys.net/10/341/2010/acp-10-341-2010-supplement.pdf the number of hourly H values higher than 200 m, N(%), determined from the observations (white bars) and from the EMEP model (blue bars) per month during the year 2001 at the Cabauw tower is presented.It should be pointed out that in this work the Ri B numbers are estimated differently from the observations and from the model.From the observations Ri B numbers are estimated using values at 2 m as the lowest level, z 1 =2 m, while Ri B estimated from the EMEP model use the first model level (z 1 ≈50 m) as the lowest level.As a consequence, considerably more cases, say ∼30%, with H >200 m are found in the observations than in the model (Fig. S1: http://www.atmos-chem-phys.net/10/341/2010/acp-10-341-2010-supplement.pdf), which is in agreement with the findings of Vogelezang and Holtslag (1996).The annual course has two maxima during spring and autumn N ∼80% in the observations and N∼70% in the model (Fig. S1: http://www.atmos-chem-phys.net/10/341/2010/acp-10-341-2010-supplement.pdf).During the winter N is expectedly smaller with N ∼60-70% from the observations and N∼30-40% from the model.During the summer N∼70-80% of cases with H >200 m is found in the www.atmos-chem-phys.net/10/341/2010/Atmos.Chem.Phys., 10, 341-364, 2010  observations and N∼50-60% from the model.Furthermore, in Fig. S2: http://www.atmos-chem-phys.net/10/341/2010/acp-10-341-2010-supplement.pdfrelation between the r and N determined from the model is shown.Obviously, N is related with r in the way that an increase in N is reflected in a decrease in r.

Significance tests
The higher level of significance for NO 2 is found at stations in Germany, Ireland, Netherlands, Norway and Sweden (not shown).Changes in r are significant over Denmark and Spain, while for SO 2− 4 there is no significant change in r with the change of the K(z) scheme in the model.other months the level of significance is satisfactory while for February and June the change in r is not significant.New parameterization schemes for K(z) and H give somewhat better results and improvements are evident although standard significance tests do not reflect it completely due to their own stated limitations in the application to this particular data.

Conclusions
Two changes of the turbulence parameterization for the EMEP model are introduced: the replacement of the Blackadar (1979) in stable andO'Brien (1970) in unstable turbulence formulations with an analytical vertical diffusion profile called Grisogono (e.g.Grisogono and Oerlemans, 2002), and a different mixing height determination, based on the bulk Richardson number formulation.The integral empirical constants, C(K) and C(z max ), are determined from the LES data in neutral and stable conditions and universally applied in the Grisogono approach for all stability conditions.The evaluation of the model performance on r and BIAS is conducted for the operational and new model setup at all available measurements from EMEP stations in the year 2001.
The main conclusions are: -The EMEP model shows a moderate improvement in r for NO 2 and SO 2 and a slight improvement for SO 2− 4 for the most of the analyzed stations.The improvements in the model performance (0.001≤D(r)≤0.12)are found at 51% of stations for NO 2 (mainly at stations in Central Europe), at 54% for SO 2 and at 55% for SO 2− 4 .However, a decreased performance (−0.09≤D(r)≤−0.001) is seen at 35% of stations for NO 2 , at 24% of stations for SO 2 in Scotland and in the shipping area, and at 22% of stations for SO 2− 4 .The annual r between the measured and modelled daily surface concentrations show slight improvements from 0.63 with the STD scheme to 0.65 with the Grisogono scheme for NO 2 , and from 0.55 to 0.57 for SO 2 .For the SO 2− 4 the correlation coefficient is around 0.61 with both schemes.
-Stations that are more affected by the local emission sources, as well as mountain stations do not show significant improvement with the change of the K(z) scheme.
On those stations the magnitude of the error was much higher than the magnitude of the variability resulting from the change of the K(z) scheme.These results indicate that a higher horizontal resolution, as well as better defined emissions, is needed in order to be able to simulate air pollution transport in complex coastal terrain under the influences of local sources.
-The new integral empirical coefficients for the K(z) in Grisogono scheme are derived from the LES data and a better balance between the momentum and sensible heat is achieved.The newly extended approach summarizes and extends the previous work (JV09) and the new integral empirical coefficients are recommended for a further model development.It is generally known that all empirical constants posses uncertainty which is affecting the accuracy of K(z) schemes.Here the accuracy of the empirical constants depends on the reliability of the LES data.Nevertheless, LES data are valuable and easily obtained data in controllable and properly idealized environment which can and should be used for model evaluation and empirical coefficients determination purposes.Further classification of the empirical coefficients according to the stability classes including convective conditions is foreseen in future work with the EMEP and EMEP4HR models.
-The Grisogono scheme is a non-local approach and it mainly depends on the position and intensity of K max .
The value of K max explicitly includes u * and H in this method, utilized from the meteorological driver and its accuracy is constrained with the NWP model performance.In air quality modelling, all K(z) schemes depend on the capabilities of used meteorological drivers as well as on model's horizontal and vertical grid resolution.Improvements in the NWP model performance would yield to appreciable differences in terms of both the magnitude and spatial distribution of pollutants which would in the end improve the air quality model performance.The Grisogono method is technically convenient since only two input variables are demanded instead of four.Therefore, the Grisogono scheme for K(z) determination is recommended for practical applications in the model yielding an improvement in overall model results.The future implementation of the integral empirical coefficients in the CBL conditions will additionally contribute to the better model performance.
-The Blackadar method, applied in the model for stable conditions, is based on the M-O theory (Monin and Obukhov, 1954).There are many studies which show that the surface-layer formulations based on the M-O theory are often not applicable in statically stable conditions (e.g.Mahrt, 1999;Pahlow et al., 2001;Poulos and Burns, 2003;Mauritsen et al., 2007;Grisogono et al., 2007).Furthermore, the Blackadar local approach generally severely overestimates the intensity of vertical turbulent transport in the LES, while it shows good results for vertical turbulent heat transport during neutral and nocturnal statically stable conditions with an enhanced turbulent mixing (K values over 30 m 2 s −1 ) in the atmosphere.It can be mentioned that such intensified local turbulence events, which the Blackadar method performed well, are a characteristic of e.g.residual layer, low level jet or clear air turbulence above the ABL (e.g.Stull, 1988), but those were not studied here in details.
-The EMEP model is able to reproduce the spatial and temporal variability of H , with r≈0.7-0.9, calculated between the model and radiosoundings during convective conditions, and r≈0.4-0.6 in statically stable conditions.It is found that the new ABL height scheme, based on the Ri B number, performs better in statically stable conditions compared to the method based on the Blackadar K(z) profiles, while the standard method has better agreement in convective conditions.The ABL height calculated with the EMEP model is generally in good agreement with the radiosounding measurements from different stations in Europe.However, due to still relatively low model resolution the model was not able to reproduce H well in complex coastal orography i.e. at Thorshaven station (Fig. 10c).At Lisbon station, which is located near the boundary of the EMEP model domain, the modelled results are dominated by weakly varying boundary conditions (Fig. 10d).It is shown on the analyzed episode of high pressure system movement across Northern Europe that accuracy of the simulated H is constrained with the NWP model performance (Fig. 8).
-The considerable number of cases with H >200 m, i.e.N, during the CBL conditions at the Cabauw tower is found, and also a negative effect of N on r values is established.The sensitivity of the Ri B scheme on the choice of the lowest layer is confirmed in this paper, showing that in the case of strong surface influenced lowest layer, a considerably more cases ∼30%, with H >200 m are found, which is in the agreement with the results of Vogelezang and Holtslag (1996).In this paper the model's ability to simulate the time evolution of the ABL and the strength of turbulence in the lowest part of the ABL is investigated and validated.Measurements at higher levels will help to identify the differences between the two ABL height schemes performances.Nevertheless, generally higher r and the similar performance of both ABL height schemes, during the warmer part of the year is in agreement with the radiosoundings results, which showed that the ABL height scheme based on the Ri B number method performs better in stable conditions than the operational one.
The evaluation study of different K(z) and ABL height schemes applied in the EMEP model provides a basis for a further model evaluation and development within the frame of the EMEP4HR project.

Figure 1 .
Figure 1.Stations used for the evaluation of the EMEP model performance.The station altitude is represented with different colours ranging from less than 300 m (blue) to higher than 3000 m (red).

Fig. 1 .
Fig. 1.Stations used for the evaluation of the EMEP model performance.The station altitude is represented with different colours ranging from less than 300 m (blue) to higher than 3000 m (red).

Figure 2 .
Figure 2. The average values and standard deviations of the new empirical coefficients for momentum a) C(K m ) and b) C((z max ) m ) and for heat c) C(K h ) and d) C((z max ) h ) calculated from the LES data.On the x-axis the number of the LES run is given, i.e. n run .

Fig. 2 .
Fig. 2. The average values and standard deviations of the new empirical coefficients for momentum (a) C(K m ) and (b) C((z max ) m ) and for heat (c) C(K h ) and (d) C((z max ) h ) calculated from the LES data.On the x-axis the number of the LES run is given, i.e. n run .

Fig. 4 .
Fig.4.Same as Fig.3.but with added K(z) profiles calculated with the Blackadar method (dashed black line), K B .

Fig. 5 .
Fig. 5.The spatial interpolation of the correlation coefficients, r, determined between the measured and modelled (a) r(NO 2 ), (b) r(SO 2 ) and (c) r(SO 2− 4 ) with the operational EMEP model, and the spatial interpolation of the differences in r, D(r), resulted due to the Grisogono scheme employed in the model for (d) D(r(NO 2 )), (e) D(r(SO 2 )) and (f) D(r(SO 2− 4 )).The available measurements from the EMEP network are used in the year 2001.

Figure 6 .
Figure 6.The annual course of: a) r, b) BIAS, c) RMSE between the measured and modelled ) ( 2 NO c with two different K(z) schemes employed in the EMEP model, i.e. the STD (red) and Grisogono scheme (blue), and d) the modelled and observed (green) monthly averages of ) ( 2 NO c in the 2001.

Fig. 6 .
Fig. 6.The annual course of: (a) r, (b) BIAS, (c) RMSE between the measured and modelled c(NO 2 ) with two different K(z) schemes employed in the EMEP model, i.e. the STD (red) and Grisogono scheme (blue), and (d) the modelled and observed (green) monthly averages of c(NO 2 ) in the 2001. 7c Fig 7. Monthly: (a) r, (b) BIAS and (c) average calculated between the ABL height ( H ), determined from the soundings ( sond H ) red bars, and H calculated from the EMEP model with the operational scheme ( std H ) grey bars and with the Ri B scheme ( new H ) blue bars for different radiosounding stations in Europe (Table 1) in January 2001 at 00 UTC.

Fig. 7 .
Fig. 7. Monthly: (a) r, (b) BIAS and (c) average calculated between the ABL height, H , determined from the soundings (H sond ) red bars, and H calculated from the EMEP model with the operational scheme (H std ) grey bars and with the Ri B scheme (H new ) blue bars for different radiosounding stations in Europe (Table 1) in January 2001 at 00:00 UTC.
) simulated in the model which occurred from 13 to 20 January 2001.The simulated lower values of H EMEP are connected with the high pressure

Fig. 11 .
Fig. 11.The monthly vertical profiles of the average hourly Ri B number calculated from the Cabauw data in from January (a) to December (l) in the year 2001.The ABL height, H , is represented with Ri B C =0.25 (thick blue line at the top of the blue area).

Fig. 12 .
Fig. 12. Monthly r between H calculated from the measurements and H calculated with the standard (H std ) -red, and the new ABL height scheme (H new ) -blue, in the EMEP model for the year 2001.

Table 1 .
The list of radio sounding stations over Europe used for validation of the ABL height, H , from the EMEP model in January and July 2001.Station name, coordinates, country, station altitude (m) and observational terms according to UTC are given.
) show a diurnal variation of H from the nighttime low H in the statically stable conditions towards the high daily H values in the convective unstable conditions.The model captures H sond daily variations and a good agreement between the H sond and H EMEP at 00:00 and 12:00 UTC is found e.g. for Meiningen r=0.91 and Madrid r=0.84 with the new ABL height scheme.Note that at Lisbon and Torshaven, H sond >>H EMEP .The modelled H EMEP values were almost constant in time and consequently the corresponding lower r and higher BIAS values are found at those stations.Note that BIAS at Lisbon is the highest among all analyzed stations.Lisbon station is located near the boundary of the EMEP model domain where the modelled results are dominated by weakly varying lateral boundary conditions.Furthermore, the model was not able to reproduce variability shown in H sond both in January and July at Torshaven station located on the Faroe Islands in the Atlantic Ocean.The Faroe Islands are situated entirely within one grid cell in the model and the model was incapable to realistically represent H in the complex coastal orography due to still relatively low model resolution.
The same procedure has been applied on r calculated between the H determined from the radiosoundings and the Cabauw data and the corresponding H values estimated with the EMEP model with the two different ABL height schemes.Although the change in r is not significant according to this test, based on the evaluation provided from the radiosounding data, the level of significance is improved for Gothenburg, Herstmonceuix, Zagreb, La Coruna and Madrid during January and for Stavanger, Copenhagen, Wroclaw, Meiningen, Vienna, Payerne and Practica di Mare in July (not shown).The change in r for Cabauw is significant during March and April; for Figure 12.Monthly r between H calculated from the measurements and H calculated with the standard ( std H ) -red, and the new ABL height scheme ( new H ) -blue, in the EMEP model for the year 2001.