The critical bulk Richardson number in urban areas: verification and application in a numerical weather prediction model

Dispersion models require hourly values of the mixing height (H) that indicates the existence and vertical extent of turbulent mixing. Urban areas, which are usually industrial areas too, have H higher than rural areas, and commonly used methods for deriving H should not be applied under the same conditions as in homogeneous conditions. The bulk Richardson number (RiB) methodwas applied to determineH over Zagreb, Croatia. Impact of urban areas on the choice of critical values of bulk Richardson number (RiBc) was explored, and different values were used for convective boundary layer (CBL) and for stable boundary layer (SBL). Aire Limitee Adaptation Dynamique development InterNational (ALADIN), a limited area numerical weather prediction (NWP) model for short-range 48-h forecasts, was used to provide one set of input parameters. Another input set comes from radio soundings. The values of H, modelled and based on compared measurements, and the correlation coefficient as well as standard deviation and bias were calculated on a large data set to determine RiBc ranges applicable in urban areas. It is shown that RiB method can be used in urban areas, and that urban RiBc should have certain limitations despite of a wide spectrum of practical values used today. Significantly increased RiBc values in SBL were determined from the NWP and soundings data, which is the consequence of increased surface roughness in the urban area. The verification of ALADIN through the determination of H was also done. Decoupling from the surface in the very SBL was detected as a consequence of the flow cease resulting in RiB becoming very large.


Introduction
Dispersion models used in air pollution meteorology have been applied on urban areas, which are characterized by anthropogenic sources of heat and pollution (e.g. Arya, 1999). In accidental situations in urban areas, modelled pollutant concentrations are required. In urban areas, mixing height (H) is modified by urban heat island, heat storage and the increased roughness due to buildings even in the absence of topographic effects. There is evidence that the urban H (H Urban ) is considerably higher for the nocturnal stable boundary layer (SBL) in comparison with the rural H (H Rural ) (Baklanov and Kuchin, 2004). These differences, H UR = H Urban − H Rural , can be up to 40% (Arya, 1999), or 45%, that is H UR ≈ 700 m (Angevine et al., 1999) especially in SBL. The greatest temperature differences between urban and rural areas are usually observed during the night when the heat from the city can produce shallow convective mixed layer, even when rural surrounding is stable (e.g. Stull, 1988).
The majority of the methods currently used to calculate H are developed for horizontally homogeneous conditions but still applied for inhomogeneous conditions; hence there is a great need for verification of their applicability for urban, nonhomogeneous conditions. We have used bulk Richardson number (Ri B ) method to validate whether it is suitable for finding H in urban area from (ALADIN) Aire Limitee Adaptation Dynamique development InterNational (Geleyn et al., 1992) NWP model, as well as from the radio soundings. The Ri B method is the standard and widely used approach to derive H from NWP models (e.g. Sørensen et al., 1996;Fay et al., 1997). It 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. One of the weak points is a considerable uncertainty in the choice of appropriate Ri Bc , and Zilitinkevich and Baklanov (2002) give a comprehensive overview of different Ri B empirical estimates. Hence, there is a wide spread of Ri Bc values used in practice. Consequently, Ri Bc can deviate, but the question is how much and under what circumstances. For example, Sørensen et al. (1996) use Ri Bc between 0.14 and 0.24, Mahrt (1981) takes the critical interval from 0.5 to 1.0, and in Fay et al. (1997) calculations of H from the German NWP model, the critical value is 0.38. It is empirically shown that Ri Bc does not have a fixed value, universally applicable in all atmospheric conditions, CBL or SBL, or for both homogeneous and non-homogeneous surfaces. Evidently, Ri Bc also depends on surface roughness, z 0 , (e.g. Zilitinkevich and Baklanov, 2002) having higher values for higher z 0 . For example, Gryning and Batchvarova (2003) find Ri Bc around 0.03 to 0.05 to perform best over sea; moreover, for the marine boundary layer Ri Bc is generally smaller than that over the land. Maryon and Best (1992) use critical gradient Richardson number (Ri c ) from 1.3 all the way up to 7.2 (which seems very high), based on radio soundings from continental Europe at Izmir in Turkey. From their study, it is difficult to come up with any final conclusion regarding Ri Bc or Ri c . An extended stable surface layer similarity theory is in Zilitinkevich and Calanca (2000) showing that turbulence may exist at almost any Ri, meaning that there is no single Ri Bc . Finally, Ri Bc may also depend on the flow history (e.g. Stull, 1988). In practical applications, we must determine the level of significant turbulence knowing that eddies do not vanish completely but get progressively smaller and less intense, hence, their influence on the mixing process becomes nearly negligible.
Assuming that Ri Bc can vary from the theoretical 0.25 value, the goal here is, upon statistical parameters obtained from a significant data set, to find the most suitable interval, i.e. the range of Ri Bc 's providing optimal H Urban values. Our intention is to find out how urban areas affect Ri Bc in CBL and SBL conditions and to investigate the possibility of using ALADIN NWP model for H calculation. We will use different nighttime and daytime Ri Bc 's applied on the radio soundings and on the model data in the period of 1 year in Zagreb, Croatia. Calculated H, determined from the model (H ALADIN ) and those determined from the soundings (H soundings ), will be analysed separately for CBL and SBL conditions. Further, a combination of Holzworth (1964Holzworth ( , 1967 method and profiles discontinuity (temperature and dew point) subjectively determined is also used to obtain mixing heights (H Holzworth ) at 12 and 00 UTC radio sounding data for the same period and location; then the H Holzworth are correlated with H ALADIN and H soundings . By determining the urban Ri Bc useful in applications, the minimum deviation in H is provided, which consequently implies a higher accuracy in dispersion modelling.

Determination of the mixing height
Both the radio soundings and NWP model data from 10 April 2003 to 10 April 2004 are used here to determine H with a commonly used Ri B method for CBL and SBL. The chosen daytime Ri Bc = 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.7, 0.9 and 1; for the nighttime Ri Bc = 0.5, 0.7, 0.9, 1, 1.5, 2 and 3. Needless to say, measurements and model deal with different but complementary approaches. While the measurements have their problems and errors (e.g. accuracy, sampling, representativity), the models are limited by their basic assumptions, resolution and parameterizations. Also, H is calculated from the radio soundings using the Ri B method and the Holzworth method. The basic idea of the Holzworth (1964Holzworth ( , 1967 method is to follow the dry adiabate from the surface with the measured/expected (maximum) temperature up to its intersection with the temperature profile from the most recent radio sounding. An addition to this procedure is a subjective determination of H from the vertical profiles (wind, temperature and humidity), used here as a refinement of the Holzworth scheme for 12 UTC data and as a basic procedure for H calculation from 00 UTC data.

Numerical weather prediction model ALADIN
Project ALADIN was proposed by Météo-France in 1990 as cooperation with the national meteorological Services of Central and Eastern Europe (Geleyn et al., 1992) in the field of NWP modelling that provides the basis for the forecasting tools of modern meteorology. ALADIN is a spectral hydrostatic limited area NWP model for short-range forecasts. Primitive prognostic equations are solved for the wind components, temperature, specific humidity and surface pressure using the two-time-level, semi-implicit, semi-Lagrangian integration scheme.
The initial and boundary conditions for the ALADIN Limited Area model for Central Europe (LACE) domain are obtained by interpolation from the analyses and forecasts of the global model ARPEGE (e.g. Marku and Fischer, 1999). During the model integration, the meteorological fields are adapted to the high resolutions of local (surface) characteristics. The spectral model with primitive equations uses two-dimensional Fourier series. Digital filter initialization method is used to filter out high-frequency buoyancy waves. Physical parameterization procedures include vertical diffusion parameterization (Louis et al., 1981) with shallow convection (Geleyn, 1987), turbulence, microphysics, soilatmosphere interaction, etc. Turbulence parameterizations differ depending on stability conditions. For the neutral surface layer (10% of the boundary layer height), turbulent fluxes are constant in the lowest few metres of the atmosphere. The cases of stable or unstable surface layer are derived from the neutral case through Monin-Obukhov's similarity theory. Equations are related to each other by numerically computed functions (Louis, 1979;Louis et al., 1981). Above the surface layer, the eddy fluxes are computed as diffusive fluxes, following the K theory. To keep continuity with the surface fluxes, the dependence of K on stability is assumed. This is a typical NWP model suite.
In the ALADIN operational version employed, there was no specific urban parameterization scheme used but the model has increased roughness over urban areas. The Croatian domain contains 127 points in the x and 109 points in the y direction (or 144 in Tellus 58A (2006), 1 x and 120 in y, with an extension zone) with an 8-km resolution in both directions. The ALADIN model solves the prognostic equations on hybrid pressure-type η coordinate (Simmons and Burridge, 1981) and the level heights are consequently not constant. Vertical resolution is not uniform and the levels are around: 17, 65, 143, 251, . . . , 1673, 1977, 2306 m, etc. This model setup is described in Ivatek-Šahdan and Tudor (2004) including an efficient dynamical adaptation to the wind field (Žagar and Rakovec, 1999).

The Richardson number method
Turbulence is the main physical process that makes mixing possible, thus characterizing the atmospheric boundary layer (ABL). Finding a parameter for practical applications that will help in a simple way to determine if the flow is laminar or turbulent is not a simple task. The flux Richardson number, Ri f , related to the production of turbulent kinetic energy (TKE), is where g is acceleration due to gravity,θ is a reference potential temperature, w θ is a vertical kinematic eddy heat flux, −u v and −v w are vertical kinematic eddy fluxes of the u and v momentum, ∂u/∂z and ∂v/∂z are the vertical gradients of the u and v wind component, respectively. For unstable flows Ri f is negative, for neutral flows it is zero and for statically stable flows, R if is positive. Richardson proposed that Ri f = 1 is a critical value, and at any Ri f < 1, static stability is insufficiently strong to prevent the mechanical generation of turbulence (e.g. Stull, 1988). For negative values of Ri f , buoyancy generation of turbulence is a dominant process. Employment of flux-gradient scheme for parameterization of turbulent fluxes yields gradient Richardson number, Ri: Substituting for the gradients in (2) with the finite differences in the coordinate frame used leads to Note that (3) is not commonly used in practice. Instead, the formulation proposed and applied by Mahrt (1981), Geleyn (1987), Troen and Mahrt (1986), Sørensen et al. (1996), Seibert et al. (2000), Batchvarova (2002, 2003), Zilitinkevich and Baklanov (2002), etc. is the bulk Richardson number: . . , 37 are the model levels. (4) Here θ 1 is potential temperature at the lowest model level z 1 andθ is the average potential temperature between levels 1 and j. The wind speed at the surface is taken as zero. The main weakness of (4) describing turbulence effects is that it might not sense the relevant, most energetic eddies and thus contribute to the model errors. Moreover, turbulence does not have to be a single valued function of only Ri. The main reason for practical use of (4) is in definition of H, which is an integral property that relates surface processes to upper processes in ABL and thus should embed non-local effects. The surface is assumed to be the main source of turbulence, with the fluxes mainly driven by the surface heat and friction. When calculating H we actually determine the height of the layer which is under direct surface influence. This is why Ri B does not have a local character.
The Richardson number formulations, Ri, and Ri B , are both used for estimations of H in practical applications. Here Ri method is also tested against Ri B method (not shown). The Ri number is very sensitive to small fluctuations in the profile, and this is avoided by using the Ri B number. Since even the Ri B method underestimates H in general (Baklanov and Kuchin, 2004), there is no reason for using the Ri method here.

Urban Zagreb area
The Zagreb area, larger than 8 × 8 km, was chosen for the comparison of the urban H calculated from the measurements and from the ALADIN NWP model. The radio sounding station Zagreb-Maksimir is in the eastern part of the city and can be considered as an urban site. Measurements in Maksimir are performed since 1956, and at first that area was basically rural with the town centre to W, and forest to NE. Due to the city growth during the last 50 yr, this area consists of factories, buildings and traffic roads nowadays.
The area is relatively flat with the mountaintop Sljeme 1 km high, 10 km to N from the site. There are around 30 buildings within a couple of hundred of metres from the site. There are buildings 4-7 m high and 30 m to W, a factory complex 20-30 m high and 200-250 m to SW. The observatory building is 70 m to S from the measuring site. Another factory is to SE, and to NE there are buildings about 40-45 m high. Park Maksimir, which is mainly forest over several km 2 , is to WNW. To N and NW there are buildings about 15-20 m high and distanced about 400-500 m. The estimated surface roughness, according to Wind Atlas Analysis and Application Program WA S P-User's Guide , for this area of a few square kilometres, is in the interval of 0.8 m ≤ z 0 ≤ 1 m.
2.3.1. Radio sounding data. Radio soundings are a common source of data for estimating H operationally. They measure all atmospheric variables of interest (pressure, humidity and temperature as well as wind speed and direction) across the full vertical profile. Unfortunately, these measurements are only taken twice a day at 00 and 12 UTC, and consequently the soundings can only be used as an overall reference. Radio sounding data are collected every 10 s in the first 5 min, every 30 s in next 15 min and every 1 min until the end of sounding path. Ascent velocity is 5 ms −1 . Crossing ABL along a skewed path within a few minutes provides a 'snapshot'-like profile that can be a limiting factor in estimating H (e.g. Parlange and Brutsaert, 1989). If the time response of the radiosonde instruments is too long, some sharp changes in the profiles can be smoothed artificially.
Large differences might occur between the depth of ABL inferred from the soundings and the depth of the turbulence especially for neutral and stable conditions. Furthermore, error in radio soundings measurements contains two components: (a) the soundings include horizontal structure in addition to the main vertical structure and (b) the sounding is almost instantaneous and therefore subject to random errors due to fluctuations caused by eddies. The contamination of the sounding by horizontal structure could be significant in an urban area where horizontal gradients may be high due to surface heterogeneity.
For calculations of H, measurements at standard and significant levels are used. Data are interpolated where missing to have the temperature, wind and height at all standard and significant levels. A simple linear interpolation method is used for the temperature, wind and height profile. For low H, there is always a large uncertainty in its determination (e.g. Grisogono et al., 1998;Mahrt, 1998) and raw data would be very useful. However, only measurements at standard and significant pressure levels are stored and archived. These measurements have typical resolutions from 100 to 1000 m and uncertainty of H is where z is the resolution with approximate values of 50 m for the first 5 min, 150 m for the next 15 min and 300 m to the end of sounding path. This suggests that the errors are smaller in the lower atmosphere, less than 150 m. Significant levels are selected from portions of the sounding where the vertical profile of temperature or relative humidity varies appreciably from a straight line. In other words, temperature and relative humidity are assumed to vary linearly with height between significant level data points. The wind direction and speed data are functions of height for 'significant' levels. Hence, wind variations are assumed to be linear between significant levels. Thus, data on standard and significant levels are acceptable for finding H with Ri B method.

Results and discussion
Using (4), two-dimensional fields of predicted Ri B (z, t) are calculated at all model levels, for every model run for this period, in the specific model grid point of Zagreb, Croatia. Time series of mixing heights, H(t), are obtained by this procedure for t = 0, 1, 2, 3, . . . , 48 h of the model output. For the same periods at the same location, H is calculated from the data provided by the radio soundings performed at 00 and 12 UTC every day. Finally, all produced data sets at the same point for the same period are correlated, and standard deviations as well as biases are calculated. This procedure yields to optimal Ri Bc for the urban area.

Verification of H using different critical bulk Richardson numbers
Our basic task here is to find the best Ri Bc for the urban area that is Ri Bc soundings when applying the bulk Richardson method on the radio sounding, and Ri Bc ALADIN , on the NWP data. In Fig. 1, daytime (a) and nighttime (b) correlation coefficients are represented. The optimal daytime Ri Bc , with the maximum correlation coefficient, for the soundings is 0.2 ≤ Ri Bc soundings ≤ 0.3, and for the model is 0.2 < Ri Bc ALADIN < 0.5. According to these results, ALADIN underestimates H soundings . However, the choice of higher Ri Bc ALADIN appears limited (see Fig. 1a). For Ri Bc ALADIN 0.7, the correlation decreases with Ri Bc soundings increasing towards 1. Nevertheless, for Ri Bc ALADIN ≥ 0.7 the correlations appear to be quite insensitive to Ri Bc soundings values. In urban CBL conditions, Ri Bc soundings does not deviate much from 0.25, the theoretical value. In SBL, Ri Bc is higher for both the NWP and sounding data (Fig. 1b). The shift to the higher Ri Bc from the day to the night is obvious. Figure 1b tells that, in average, the more stratified nighttime turbulence, occurring at relatively smaller scales than its daytime counterpart, requires larger Ri Bc . This is in accord with, e.g. Stull (1988) and Zilitinkevich and Calanca (2000) and also implying a hysteresis effect in various critical Ri numbers. The optimal nighttime critical values are 1.5 < Ri Bc soundings ≤ 2, and 1 ≤ Ri Bc ALADIN ≤ 1.5. The scatter plots with Ri Bc soundings and Ri Bc ALADIN giving the maximum correlation coefficient are shown in Fig. 2. For CBL, the maximum correlation coefficient r CBL = 0.84 is obtained for (Ri Bc soundings , Ri Bc ALADIN ) day = (0.2, 0.4), and for the SBL it is r SBL = 0.86 for (Ri Bc soundings , Ri Bc ALADIN ) night = (1.5, 1.5).
To find the best combination of the Ri Bc 's that should be used in practice, we take into account standard deviations and biases. It is expected that the biases and deviations are small for the (Ri Bc soundings , Ri Bc ALADIN ) combination that has the highest cor-   Fig. 3. Figure 3a confirms that the minimum standard deviation of 271.5 m is found for (Ri Bc soundings , Ri Bc ALADIN ) day = (0.2, 0.4). The nighttime standard deviation for (Ri Bc soundings , Ri Bc ALADIN ) night = (1.5, 1.5) is 227.5 m (Fig. 3b). The minimum deviation is found for 0.5 ≤ Ri Bc < 2.5, for both the sounding and ALADIN data.
The absolute bias, |BIAS| = |H sounding −H ALADIN |, is shown for again daytime (Fig. 4a) and night time (Fig. 4b) conditions. In CBL, the minimum bias is for 0.2 ≤ Ri Bc soundings < 0.5, and 0.3 < Ri Bc ALADIN < 0.6. Furthermore, the bias for SBL is smallest for 0.5 < Ri Bc soundings < 2 and 1 < Ri Bc ALADIN < 2. These results confirm good selection of Ri Bc ALADIN and Ri Bc soundings shown in Fig 2. As expected, for the maximum correlation coefficient the corresponding standard deviation and bias minimize. Based on the results shown in Figs. 1 through 4, the decision for the optimal Ri Bc can be made.

The Holzworth method
Another method is used to objectively analyse the results and to find out the accuracy of the Ri Bc selection in urban areas. Here H from the sounding is determined with the Holzworth method. This H Holzworth is correlated with H ALADIN and H sounding calculated with the bulk Richardson method for different Ri Bc values. The goal is to find the best Ri Bc ALADIN and Ri Bc soundings for CBL and SBL conditions based on the maximum correlation coefficient and the minimum standard deviation and bias determined between H ALADIN and H sounding or H Holzworth . The results are in Fig. 5 and they are compared with the corresponding results in Figs. 1 through 4.
The maximum correlation in CBL for soundings is achieved for 0.1 ≤ Ri Bc soundings ≤ 0.3. These results from Fig. 5 are compared with those in Fig. 1a. The optimal Ri Bc soundings is almost the same. This confirms that in CBL conditions Ri Bc soundings should not deviate much from the theoretical value 0.25, even in urban areas. It means that when buoyant processes dominate in generation of turbulence, other impacts on Ri Bc are less relevant for H calculations. In SBL conditions, the best results are found for 0.7 ≤ Ri Bc sounding ≤ 2 and the shift to the higher values is confirmed. The optimal value is Ri Bc sounding = 1.5, same as before, and the corresponding bias is only 13.2 m.
The Ri Bc ALADIN results confirm two important earlier conclusions for CBL. These are a need for higher Ri Bc in the NWP model, here 0.4 ≤ Ri Bc ALADIN ≤ 0.7 is found, and significant decrease of the correlation coefficient for Ri Bc ALADIN > 0.7. Note that the correlation is smaller and less sensitive for the NWP data in Fig. 5 than in Fig. 1a, which is a consequence of the different methods (bulk Richardson and Holzworth method) employed for the determination of H. This is also confirmed with the correspondingly higher standard deviation and bias. Also note that for CBL, the NWP model standard deviation is around 500 m and less sensitive for different Ri Bc values. Furthermore, bias values in CBL for model and for soundings give only tendency to higher values for higher Ri Bc . In SBL, the correlation coefficient and standard deviation give 0.7 ≤ Ri Bc ALADIN ≤ 2, which is the same as for the soundings and in agreement with previous results. Nevertheless, bias shows a tendency to higher Ri Bc ALADIN values having its minimum for 1.5 ≤ Ri Bc ALADIN ≤ 2.5.

Vertical variations of Ri B and frictional decoupling in the NWP data
Here vertical daily and seasonal variations of Ri B , ∂ Ri B ∂z are analyzed; occurrence of frictional decoupling (FD) that may lead to numerical instabilities, is found. Results for three winter days (Fig. 6) and for three summer days (Fig. 7) represent typical daily vertical evolution of Ri B starting at 00 UTC. Two critical values are chosen, Ri Bc = 0.3 and 1, marked on the figs with thick black curves, illustrating H and showing the model sensitivity to the seasonal variation acting on H. At the switch into a new day, the curves in Figs. 6 and 7 are discontinuous due to the model re-initialization.
In Fig. 6, we see a realistic daily evolution of H. The difference between H determined with Ri Bc = 0.3 and 1 is around 100 m for SBL, while for CBL this difference can be more than 1000 m. Once again, this proves the importance of the correct selection of the Ri Bc . The synoptic situation related to 4 and 5 February was characterized by long stable periods with fog and low stratified clouds due to high-pressure field over Croatia and warmer NW airflow from the Atlantic Ocean at higher altitudes. However, the stable synoptic situation ended on 6 February with a Genoa cyclone approaching and intensifying SW wind. The ∂ Ri B ∂z is also different for 6 February than on other two days.
Summer days in Fig. 7 show higher H, as expected. An elevated cyclone formed west from Zagreb determining the synoptic situation on 13 August 2002. There was advection of moist, unstable and relatively warm air from southwest. On 14 August, the synoptic situation stabilized. During the next day, 15 August, the weather was sunny and warm. The model gives a very intensive convection later in the evening and night Tellus 58A (2006), 1 which can be a prediction error. Exception of these data would slightly increase the correlation; here results contain model errors too.
Note, a nearly collapsed SBL around midnight between 4 and 5 February in Fig. 6, and for summer days around midnight in Fig. 7. This emphasizes one of the unsolved turbulence parameterization problems when the flow ceases, U ≈ 0 and Ri → ∞ (e.g. Mahrt, 1998;Zilitinkevich and Calanca, 2000;Grisogono and Oerlemans, 2001), producing FD in almost any NWP model and ALADIN as well. The FD is connected with an increase of stability when the flow laminarizes near the surface. It occurs when shear/friction fails as the dominant generator of turbulent fluxes.

Conclusions
It is shown that it is necessary to provide adjustments of the bulk Richardson method via the critical bulk Richardson number, Ri Bc . These adjustments are needed in urban areas to calculate H based on measurements and from NWP models.  As it may have been expected, when applying the bulk Richardson method on the radio sounding data in CBL, Ri Bc soundings does not deviate much from the theoretical value even in the urban area. However, in SBL, where the dominant eddies are relatively small, Ri Bc soundings should be higher, namely here 1.5 Ri Bc soundings ≤ 2. Such larger Ri Bc indicates the urban influence; increased surface roughness enhances Ri Bc soundings . This is also confirmed with the results calculated with the Holzworth method in Fig. 5.
The results from ALADIN data in CBL show a little more spread towards larger R iBc . Based on H derived with the bulk Richardson method as well as with the Holzworth method, we conclude that in CBL Ri Bc ALADIN ≤ 0.7 when applied on NWP data. The optimal CBL Ri Bc ALADIN = 0.4 corresponds to the correlation coefficient r CBL = 0.84. Higher CBL Ri Bc ALADIN tells us more about NWP model features and indicates that urban parameterizations ought to be included in the model. This should not be taken straightforwardly, and further studies and improvements in the model are needed. In SBL, 1 ≤ Ri Bc ALADIN ≤ 1.5 is found optimal when H ALADIN is correlated with H soundings determined with the bulk Richardson method and 0.7 ≤ Ri Bc ALADIN ≤ 2 when H soundings is determined from discontinuities in the vertical profiles.
The maximum correlation coefficient in SBL is r SBL = 0.86 for the (Ri Bc soundings , Ri Bc ALADIN ) night = (1.5, 1.5). It is shown that Ri Bc increases with an increase of stability. Nevertheless, in very stable conditions, the model FD from the surface is detected.
The persisting FD may result in numerical instabilities, and it must be excluded from the model. The FD effect indicates that surface parameters must be included when calculating H from an NWP model, especially in SBL.
It should be pointed out that the correlations determined here contain the model and sounding errors too, and the results may not be interpreted only as methodology insufficiencies. Also, observed ABL depth is not based on actual turbulence measurements and is itself a rough approximation, especially at night.
Based on more than 300 cases of H derived from the measurements, it is found that Ri Bc > 0.25 in urban SBL; here optimal value is 1.5. Hence, larger Richardson number is required with greater roughness of the urban area, meaning that for a given large value of the stratification, turbulence is more likely generated over a rough surface than a smoother surface. This also implies that a simple parameter, e.g., various forms of Ri based on similarity theory, seems insufficient for determining H in urban SBL, which ought to be a multi-valued function of a couple additional parameters (e.g. Vogelezang and Holtslag, 1996;Gryning and Batchvarova, 2002;.

Acknowledgments
The authors are thankful to our colleague Vesna D uričić who diligently determined mixing heights from Zagreb radio soundings with the Holzworth method. The authors are grateful to Sven-Erik Gryning on his valuable advices and given important Tellus 58A (2006), 1 aspects concerning urban mixing height. This work was performed at the Hydrological and Meteorological Service of Croatia and supported by Croatian Ministry of Science, Education and Sport, under the projects number 0004001; B. G. thanks to the same source, project 0119330. Three anonymous reviewers gave many valuable comments. We appreciate the reviewers' constructive criticisms, which helped to improve the paper. Lukša Kraljević, KornelijaŠ.Čanić, Stjepan I.Šahdan, Martina Tudor and Zvonimir Jeričević are thanked for technical support.