A study on parametrization of orography-related momentum fluxes in a synoptic-scale NWP model

In a synoptic-scale atmospheric model, winds and subgrid-scale momentum fluxes are influenced by the properties of the underlying surface. Parametrization schemes have been developed to describe the interactions of the air flow with the rough surface and small- and mesoscale orography. The schemes interact with each other and with the model’s dynamical processes. Newparametrizations were introduced to the High Resolution Limited Area Model (HIRLAM), to account for the subgrid-scale orography effects. Needed orography parameters were derived from a high-resolution digital elevation data set. A detailed analysis of the momentum fluxes and kinetic energy budget helped to understand the interactions between parametrizations. Mean and root mean square error, averaged over all observations and forecasts, did not reveal significant differences between the updated and reference model results. However, more detailed diagnosis of the forecast-observation differences allowed to show that the new parametrizations of mesoscale and small-scale orography lead to more realistic low-level wind distribution.


Introduction
Tendencies of the horizontal wind v (x, y, z) given by a numerical weather prediction (NWP) model consist of contributions from explicitly resolved dynamics and physical parametrizations of subgrid-scale processes: where the indexes d and p refer to dynamics and parametrized physics. In a synoptic-scale model, the dynamical tendency contains terms describing advection, pressure gradient and Coriolis forces. In addition, it is influenced by the accuracy of the numerical schemes, smoothing, boundary relaxation, interpolations, etc. The parametrized tendency is due to the divergence of the stress tensor τ ij , which defines the sum of surface forces upon a fluid element. This tendency is conveniently approximated by the vertical divergence of subgrid-scale vertical momentum fluxes, where the three-dimensional vector τ (x, y, z) is the sum of momentum fluxes due to n different physical processes, denoted here by the index j, ρ is the air density, w is the vertical velocity, an overline denotes gridbox average and a prime subgrid-scale deviation.
Corresponding author. e-mail: laura.rontu@fmi.fi A classical parametrized process is the turbulence over a rough surface. Here, the properties of the surface (aerodynamic roughness) and the flow (thermal stratification, velocity) are accounted for, usually separately for the surface (constant flux) layer between the surface and the lowest model level, and for the layers higher in the planetary boundary layer (PBL) and free atmosphere. Turbulence parametrizations relate the subgrid-scale momentum fluxes to the vertical gradients of the resolved wind components.
In addition, parametrizations of the momentum fluxes due to several other processes including the effects of unresolved smallscale orography (SSO; "orographic stress"), wave drag due to the breaking orographic or non-orographic buoyancy (gravity) waves, drag due to blocking of the low-level flow by mesoscale mountains or effects of local circulations like sea breeze have been developed.
The concept of effective or orographic roughness length was developed by Fiedler and Panofsky (1972) and Mason (1985). Experimental and theoretical results of Grant and Mason (1990), Wood and Mason (1993), among others, supported the idea that over small-scale hills, the area-averaged velocity profile varies logarithmically with height in a similar way as in the small-scale turbulence. Later, the orographic roughness length has been used in the NWP models as a convenient tuning parameter, even to represent processes of larger scale than it was developed for, see for example, the model comparison study by Georgelin et al. (2000). An alternative parametrization of the SSO effects was proposed by Wood et al. (2001) and Brown and Wood (2001), see also Wilson (2003). Here, a three-dimensional momentum sink is formulated and allowed to influence directly the momentum equations of a NWP model. A practical application of this method to the ECMWF model was suggested by Beljaars et al. (2004).
The need to parametrize the wave drag in NWP models was understood already by Lilly (1972) and applied first by Boer et al. (1984) and Palmer et al. (1986). The first generation parametrizations were based on linear two-dimensional theory of mountain waves developed since 1930s (for a review see Smith, 1979). More recent formulations, for example, Lott and Miller (1997); Gregory et al. (1998); Scinocca and McFarlane (2000) aim at handling also the non-linear processes of wave generation, breaking and reflection. The blocked-flow drag parametrization by Lott and Miller (1997) uses the approach of an early laboratory study by Hunt and Snyder (1980) and also takes into account directional and asymmetry properties of the unresolved terrain (Baines and Palmer, 1990, see also Baines, 1998). A newly formulated mesoscale orography (MSO) parametrization scheme for the UKMO model (Webster et al., 2003) represents the blocked flow and wave drag by an empirical fit to the numerical simulation results of Olafsson and Bougeault (1997), where friction and rotation effects were also accounted for. An advanced orographic drag parametrization scheme was suggested by Kim and Doyle (2005). It enhances the parametrization of Kim and Arakawa (1995) to take into account the orographic anisotropy and flow blocking. Recent studies by Nappo et al. (2004) and Smith et al. (2005) are devoted to the parametrization of wave stress in the planetary boundary layer.
Each of the different parametrization schemes has been developed based on the knowledge about a particular physical process. However, for historical reasons, new parametrizations have been developed from and combined with the existing ones. Thus, the effective roughness length aimed at parametrization of the unresolved orography effects is combined with the aerodynamic roughness parameter ("vegetation roughness") derived from the land cover characteristics such as vegetation. The parametrized drag due to the flow blocking is formulated as a modification of the wave drag vector. On the other hand, a unified approach to handle wave and turbulent drag due to different scales of orography is still lacking.
Although all the parametrized processes are subgrid scale relative to a NWP model of a given resolution, the orography scales involved in the different subgrid-scale processes are different. In the parametrization scheme of MSO effects developed by Scinocca and McFarlane (2000), special care was taken to separate the orography scales treated by the parametrization and resolved by the model. A similar approach is used in the derivation of the orography-related parameters for HIRLAM  and for the recent version of UKMO model (Webster et al., 2003). A spectral approach is applied in the new parametrization scheme of SSO effects for the ECMWF model (Beljaars et al., 2004).
The purpose of a numerical weather prediction model is to simulate the state of the atmosphere as a whole, including all possible interactions between the processes described by the separate blocks of the model. However, the relations, interactions and scales relevant for the different processes are not always understood and handled systematically. As noted by Kim and Doyle (2005), the distinction among various drag mechanisms near the surface has been rather vague in the literature. Experiments with the parametrization scheme for MSO effects of the Météo France ARPEGE model (Geleyn et al., 1994), also adapted for HIRLAM, showed that the turbulence and MSO parametrization schemes are strongly coupled and compensate each other in model simulations . In a global model, similar behaviour was revealed by Kim and Hogan (2004).
The aim of the present study is to analyse systematically the influence of the orography-related parametrizations of momentum fluxes on the simulated low-level wind and understand the mutual interactions of the parametrized processes. A limited-area atmospheric research and NWP model HIRLAM (Undén et al., 2002) is used as the basic tool and reference model of the study. New parametrization schemes for the MSO and SSO effects are introduced into HIRLAM. The derivation of orography-related variables of HIRLAM, based on digital elevation data, is renewed. We will concentrate on lower troposphere momentum fluxes in a synoptic-scale version of the model. The rest of the paper is organised as follows. Section 2 describes the parametrization schemes modified and used in this study. Section 3 presents the derivation of orography-related parameters. Description of the numerical experiments, their results and discussion follow in Sections 4 and 5. Conclusions and outlook finish the study in Section 6. Table 1 lists the parametrization schemes applied and developed in this study. The "reference" refers to the present HIRLAM and "modified" to the suggested parametrizations. A description of the modified schemes is given below. For convenience, a short summary of the relevant reference schemes of HIRLAM is also included. The typical scales of the underlying orography related to the schemes are shown in Table 2 and discussed in more detail in Section 3.

Turbulent fluxes over flat surface
Momentum flux over a (flat) rough surface is parametrized in HIRLAM in a conventional way (Undén et al., 2002). In the surface layer between the lowest model level and the surface it is defined by the lowest-level wind and a drag coefficient: where v s = 0 is the surface wind and v nlev denotes wind vector at the lowest model level. (Note that positive surface stress Tellus 58A (2006), 1 is associated with downward momentum flux.) The indexes s and t refer to surface and turbulence, respectively. The functions used for computation of the aerodynamic resistance in the surface layer are based on Louis (1979) and Louis et al. (1982). The drag coefficient C u depends on surface roughness, stability and wind shear (Undén et al., 2002). In the reference HIRLAM the momentum roughness consists of both vegetation (z o,veg ) and orographic (z o,oro ) components (see Section 2.2). In the present modified experiments, only vegetation roughness (z o,veg ) is used for calculation of the turbulent momentum fluxes. The resulting turbulent stress vector τ ts is turned clockwise (in the Northern Hemisphere) with respect to the surface layer wind v nlev in the (stable) boundary layer (Nielsen and Sass, 2004). The surface turbulent momentum flux given by eq. (3) defines lower boundary conditions for a TKE − l turbulence scheme based on Cuxart et al. (2000), which parametrizes the turbulent fluxes above the surface layer using the classical eddy diffusivity formula: The coefficient K is computed based on the turbulent kinetic energy TKE and a diagnostic length scale l (Lenderlink and Holtslag, 2004). TKE values are obtained by solving a prognostic equation, for details see Undén et al. (2002).

Drag due to unresolved small-scale orography
In the reference HIRLAM, an orographic roughness length z o,oro is calculated based on the subgrid-scale orography variance and the number of peaks within a grid-square (Undén et al., 2002). Empirical scaling parameters are used to enhance the value of this parameter over steep orography . 1 With this approach, all obstacles with a scale smaller than the size of the model's grid-square contribute to the orographic roughness. For the calculation of the surface layer momentum fluxes, z o,oro is combined with the vegetation roughness, z o,veg , derived from the properties of land cover.
In the present study, a parametrization of SSO effects based on Wood et al. (2001) and Brown and Wood (2001), see also Wilson (2003), is implemented in HIRLAM and used instead of the orographic roughness. Here, the orographic stress τ o due to the small-scale orographic features is estimated by the expression where C o is an orographic drag coefficient, s t denotes the mean maximum small-scale slope (tangent) over the grid-square, l o is a decay scale in vertical and z is the height above surface. According to eq. (5), the surface orographic stress τ os = C o τts ρs s 2 t is parallel to the turbulent stress τ ts , which is determined by the wind and stability in model's surface layer (see Section 2.1). The vertical decay of the orographic stress is determined by a simple exponential function. In eq. (5) several simplifications to the formulae proposed by Wood et al. (2001) were made. The direction differences of slopes and the effects of anisotropy inside the grid-square are ignored, because only the smallest scale (presumably isotropic) features are allowed to influence. The parameter C o was given a constant value (C o = 2000) instead of taking into account its possible dependency on wind shear, etc. The value of l o was set equal to twice the SSO standard deviation σ t in a grid-square.
Derivation of the variables s t and σ t from digital elevation data is discussed in Section 3.2.
The suggested formulation differs from that by Beljaars et al. (2004), where a spectral approach to the description of subgridscale orography height is applied. With consequent simplifications, the authors end up in an approximation for the tendency of the horizontal wind v(x, y, z) (their eq. 16 rewritten in a form analogous to eqs. 2 and 5): where C 1 is an orographic drag coefficient, f 1 (z) is an empirical function of height derived from the analysis of the spectral properties of SSO, σ 2 flt is a filtered variance of orography height (representing horizontal scales between 5 and 20 km, close to σ m defined below in Section 3.3). Here, the wind tendency at any level depends on the wind at the same level.

Mesoscale orography effects
A parametrization scheme for MSO effects, developed for the Météo France ARPEGE model, is being implemented into the reference HIRLAM system. The scheme handles vertically propagating buoyancy (gravity) waves and blocking of the low-level flow due to mesoscale mountain systems. A short summary is given here, for details see Rontu et al. (2002).
The wave drag is estimated by a formula based on the linear two-dimensional theory, where the index s refers to mean near-surface values, K g is a tuning parameter depending on the model resolution (here K g = 3.5 · 10 −6 m −1 ), N is the buoyancy (Brunt-Väisalä) frequency, N 2 = g θ ∂θ ∂z , h m is the subgrid-scale mountain height based on the standard deviation of MSO and v f s is a (fictive) surface wind representing the layer between surface and h m and parallel to the stress vector τ ws .
As long as there is no wave dissipation, the wave momentum flux is constant with height. The momentum sink is realized when the waves break. The parametrization of wave breaking processes follows Lindzen's saturation theory (Lindzen, 1981). In addition, (non-linear) wave reflection from a breaking level is taken into account. Wave breaking and reflection modify the surface value τ ws and the profile of the wave drag τ w (z).
Low-level flow blocking is assumed if a non-dimensional mountain height G (also called inverse Froude number based on mountain height, see Baines (1998) for clarification of the concepts) depending on stability, mountain height and upstream wind exceeds a critical value. G is defined conveniently as where U p is the velocity of the upstream wind component perpendicular to the ridge below h m . The blocked flow stress τ m at each (low troposphere) model level is calculated according to Lott and Miller (1997). Finally, it is combined with the wave drag vector τ w .

Derivation of orography-related parameters from digital elevation data
The basic variables needed by the parametrization schemes and used for definition of the average terrain height of HIRLAM grid-squares are derived from a digital elevation data set. The reference HIRLAM system uses the global set from U.S. Geological Survey 30 arcsecond data (GTOPO30, USGS (1998)) in its original latitude-longitude representation.
In the present study, calculation of the orography-related parameters is renewed and based on the 1 × 1 km resolution Hydro1K elevation derivative database (USGS, 2003). These data are derived from the GTOPO30 data but represented in an equal area azimuthal Lambert projection. Because of the areaconserving properties of this projection, the data set is well suited for further use in the rotated latitude-longitude coordinate system of HIRLAM (Undén et al., 2002). The Hydro1K data, prepared especially for the hydrological use, cover the whole globe excluding Greenland and Antarktis.
Only the terrain elevation of Hydro1K data is used here. From it, the variables listed in Table 2 are calculated for each HIRLAM grid-square of a given resolution. The averages are related to the orography scales indicated in the rightmost column of the table. The choice of 3 km for the limit between small-scale (turbulent) and mesoscale (wave related) processes is somewhat arbitraryin reality features of the flow such as stability and flow velocity should influence the definition. The limiting scale between parametrized mesoscale processes and those explicitly resolved by the model's dynamics is chosen to be two gridlengths (2 x). Vegetation roughness (whose values in the reference HIRLAM are derived from USGS (1997)) is included into Table 2 for completeness sake.

Mean elevation
Mean elevation (H 2 x ) is used for calculation of the surface pressure, on which the vertical hybrid coordinate of HIRLAM is based. In the present study, it is calculated through double smoothing over a HIRLAM grid-square. First, at each point in the 1-km-resolution Lambert grid, nine-point averages (h 9 ) of the source Hydro1K elevation (h) data are calculated. These averages now represent an area of 3 × 3 km. They are further averaged, in the rotated HIRLAM coordinate system, over 2 x-resolution squares centred at the middle of HIRLAM grid-squares. The rationale behind these operations is to avoid artificial small-and mesoscale features, which could be sources of numerical noise for the model integrations. As shown for example, by Davies and Brown (2001), atmospheric models are not able to explicitly represent underlying orography features smaller than about 4-6 gridlengths.
Smoothing is often achieved through scale-dependent filtering of the mean elevation data. The reference HIRLAM system uses a Raymond filter (Raymond, 1988) for this purpose. Comparison of the reference HIRLAM mean elevation with the one calculated in this study did not reveal significant differences, so the new one is used in all experiments reported in Section 4.

Small-scale orography
The basic parameters for the new orographic turbulence parametrization are the small-scale slope s t and standard deviation σ t (see eq. 5 and related definitions). These parameters are calculated based on the difference between the original (h) and mean (h 9 ) elevation values. The deviation at the point i, contains only the smallest-scale features. A maximum slope s max is calculated at each Hydro1K point using the difference of h i between the point i and its eight neighbours. The maximum slope values are averaged over each HIRLAM grid-square to get the averaged maximum slope s t . Also the small-scale variance σ t over a HIRLAM grid-square is derived from the h i values. Distribution of the reference HIRLAM orographic roughness and the new s t and σ t in two different resolutions over Iceland is illustrated in Fig. 1. Maximum values and area averages of the   Table 3. General features of the small-scale parameters are somewhat similar. As expected, the new small-scale variables s t and σ t are only mildly scale dependent, as the smallest-scale variations should not be sensitive to grid-size. Scale dependence may only arise because the averages are taken over grid-squares covering areas of different surface types. In addition, the influence of the maximum slopes found in a grid-square is distributed smoothly over larger areas when the resolution grows coarser. Note also

Mesoscale orography
MSO variables should represent the orography scales between those resolved by the model dynamics and those included in the parametrization of the SSO effects. Here, their values are derived from the difference of the averaged values h 9 − H 2 x . They are calculated for each grid-point over a square of 2 x × 2 x. After the averaging operations, no further filtering of the MSO variables is applied. The basic variable for the ARPEGE-HIRLAM MSO parametrization scheme is the standard deviation of MSO σ m (Fig. 1, Table 3). As expected, its values depend clearly on the model's resolution. The two other MSO parameters are the anisotropy α and the angle θ between model's x-axis and the principal axis (direction of maximum gradient) of the MSO. All three a) b) Fig 3. Surface stress components | τ ts | (solid), | τ os | (dotted) and | τ ms | (dashed) with the resolutions 11 km (circles on the lines) and 33 km (no markers), for the experiments O11 and O33 (a) and for the experiments B11 and B33 (b), averaged over 31 00UTC+48h forecasts along the rotated latitude 6.8 N in Southern Iceland, unit Pa. To retain only significant features, the time-averaged model output fields were spectrally smoothed before drawing the cross-section values.
variables are calculated according to the formulae explained in Rontu et al. (2002) but using the differences h 9 − H 2 x instead of the filtered source elevation data.

Definition of the experiments
To test the proposed modifications, a series of 1-month simulations in January, 2000 was run using the newest available HIRLAM version (essentially equal to 6.3.5). The two nested experiment domains with horizontal resolutions of 33 km and 11 km are shown in Fig. 2. A model configuration with 40 levels in vertical between the surface and the isobaric level 10 hPa was used. Full data assimilation with an interval of 6 h was performed in all experiments. Long forecasts with a length till +48 h were started once a day at 00UTC. Lateral boundary values were interpolated from archived HIRLAM analyses with a horizontal resolution of 33 km. During January 2000, a north-westerly flow over Northern Atlantic towards Scandinavia was prevailing. Over Iceland, the mean flow was westerly, with an orographic anticyclone forming over the high-elevation glacier Vatnajökull. During 25-27 January there was significant wave activity over Scandinavian mountains, indicated also by the observed stratospheric wave clouds (Dörnbrack et al., 2002).
The main numerical experiments are summarized in Table 4. MSO refers to the mesoscale (Section 2.3) and SSO to smallscale (Section 2.2) orography parametrizations. ECMWF SSO parametrizations are based on eq. (6), where σ m from the present 11-km resolution experiments is taken to represent σ flt and coefficients are taken from Beljaars et al. (2004), while HIRLAM SSO parametrizations use eqs. (2) and (5). The experiments are referred to as R33 or O11, etc., to indicate both the experiment name and horizontal resolution. Note that the new mean orography (see Section 3.1) was used in all experiments.

Diagnostic variables and tools
Averaging eq. (1) over long enough period of time allows to estimate the response of model dynamics to the forcing due to the physical processes. 2 The three-dimensional wind tendency (eqs. 2, 6) due to a physical process can be integrated vertically to obtain the corresponding surface stress, assuming that there are no vertical momentum fluxes at the top of the model atmosphere. Diagnostics collected at each time step of a HIRLAM experiment include accumulated tendencies of model variables (such as wind components or derived from these resolved-scale kinetic energy) plus independently accumulated surface fluxes. These allow to understand the role of different parametrized processes and the response of model dynamics.
The resolved-scale kinetic energy can be used for the analysis of the results of the numerical experiments. Its true sources and sinks represent the work done by frictional forces, that is, the parametrized turbulence, MSO and SSO effects. It can be estimated based on the values of wind components and their accumulated tendencies: where k = ρ 2 (u 2 + v 2 ) denotes the (diagnostic) resolved-scale kinetic energy and the index j denotes any of the different processes so that j = d, t, o, m refer to dynamics, turbulence, SSO or MSO parametrizations, respectively. The total tendency, due to the model dynamics and all parametrizations, is obtained from the difference of two model states. Averaged over long enough time it is close to zero so that the dynamical and parametrized tendencies balance each other.
Analysis of the results is performed over the common area of the experiments (see Fig. 2). For the comparison, results of the coarse resolution experiments are interpolated to the fine resolution grid. Special attention is paid to the behaviour of the parametrizations over the complex orography of Iceland. In addition, the standard station verification provides bias and root mean square errors of the forecasts against observations.

Surface stress
A west-east cross-section over Southern Iceland (along the line shown in Fig. 1) illustrates a typical behaviour of surface momentum fluxes or surface stress components due to turbulence, SSO and MSO (Fig. 3). Results of the experiments O11/O33 and B11/B33 are shown. Several features can be noted: (1) the effect of the MSO parametrizations decreases dramatically with improving resolution, (2) turbulent and SSO stress increase slightly with increasing resolution, (3) the role of oro- The different behaviour of the total surface stress in the experiments is illustrated by Fig. 5 for the whole Iceland in the case of 33 km resolution. The experiment B33 produces clearly larger stress than R33 and O33, the difference increasing with the increasing total stress values. The difference is less pronounced over the whole European experiment domain and practically non-existent over flat land areas (not shown). The differences are also smaller in the 11 km experiments.
Finally, Table 5 summarizes the comparison of averaged over whole Iceland parametrized surface stress due to the different processes and systematic differences between experiments of different resolutions. The numbers confirm the conclusions based on Figs. 3-5. Interactions between the parametrization schemes are further discussed in the next section.

Kinetic energy budget
Vertical profiles of the area-averaged kinetic energy tendencies are compared in Figs. 6 and 7. In the reference HIRLAM, turbulence-based on the effective roughness-is the only parametrized process responsible for the subgrid-scale change of the resolved-scale kinetic energy (experiment R33 in Fig. 6, upper panel). Its influence and the response of model's dynamics balance each other in the whole planetary boundary layer (PBL), in this case till the model level 30, that is, approximately to the height of 1.5 km. In the experiment O33 the effects of MSO and SSO parametrizations are smaller than that of the parametrized turbulence and seen in a shallower layer (Fig. 6, upper panel). As soon as these parametrizations are included, compensation between the schemes shows up (see the difference in Fig. 6,  Fig 6. Resolved-scale kinetic energy tendencies due to the effects of parametrized turbulence (turb, long-dashed), mesoscale orography (mso, dotted) and small-scale orography (sso, long-short-dashed), averaged over the area of Iceland and based on the accumulated tendencies of 31 00UTC+48h forecasts during January 2000. The response of model dynamics, equal in magnitude and opposite in sign to the sum of parametrized processes, is also shown (dyn, solid). lower left), as noted earlier by Rontu et al. (2002). The dynamic response due to the total stress divergence changes only a little.
Comparison of the different resolution experiments O11 and O33 (Fig. 6, lower right) shows, with improving resolution, a decrease of the MSO effect and a corresponding increase of the turbulence and SSO effects on kinetic energy. In this case the dynamic response also changes, that is, when the resolution grows higher, the predicted kinetic energy tendency due to the parametrized processes tends to increase at the lowest levels and decrease in the upper part of the PBL, that is, the effects concentrate on a shallower layer.
Experiments 'B' behave differently from the experiments 'O' both in the high-and the low-resolution cases (Fig. 7, upper panel for the 33-km resolution, comparison between O11 and B11 not shown). In B11 and B33 the SSO effects dominate the mesoscale effects and tend to retard the flow in a deep layer. Close to the surface the compensating effect of the non-orographic turbulence increases in 'B' (even creating resolved-scale kinetic energy!) compared to 'O'. Compared to 'O' (Fig. 7, lower left) and 'R' (not shown), the experiments 'B' seem to create the strongest dynamic response (retardation) in the upper part of the PBL. The difference between experiments B11 and B33 is mainly due to the decrease of the MSO effect as in the case of the experiment O11 compared to O33.
MSO, SSO and turbulence parametrization schemes are coupled with each other mainly through the low-level wind velocity (see Section 5.3 for discussion of the low-level wind), which influences all the surface stress components. The surface layer stability is expected to play a minor role. In the HIRLAM system, each of the parametrization schemes independently creates a wind tendency ∂ v ∂t using wind vector values at the previous time step. Before starting the next time step, these contributions are summed. They are then added to the tendency due to the model dynamics, and the values of the wind components updated accordingly.
In the model, a (non-orographic) turbulence parametrization scheme tends to smooth subgrid-scale vertical gradients independently of their origin or scale. The orography-related parametrizations and model dynamics, each within their own horizontal and vertical scales, may create such gradients. This might offer a possible mechanism of compensation between the schemes, shown in Figs. 6 and 7.

Surface pressure and low-level winds
In HIRLAM, a long-lasting feature seen in comparisons of the predicted with observed 10-m level wind (V 10m ), has been a systematic overestimation of wind speed by 0.5-1 ms −1 . In the present study, standard station verification (bias and root mean square error) against synoptic observations over European or Scandinavian land stations did not reveal significant differences between the experiments O33 and R33 during January 2000 (not shown). In the experiment B33, the typical for the reference HIRLAM positive bias of V 10m is slightly improved but a small positive mean sea-level pressure (pmsl) bias of about 0.2 hPa appears. Small changes of the 2-m temperature and humidity bias also appear. Verification of all 11 km experiments over Scandinavia shows negligible differences. Thus, the possible effects of orography-related parametrizations seem to be mostly local. In the following, verification results over Iceland are analysed more closely.
Here, instead of the familiar positive bias of V 10m , there is a negative one, of about −1 . . . −1.5 ms −1 (depending on the forecast length) in 'R' and 'O', and about −2 ms −1 in 'B', with a corresponding positive pmsl bias-the smallest in 'R' -of the order of 0.2-0.8 hPa (comparison not shown). Thus, the nearsurface wind seems to be retarded too much in all experiments, and the unmodified reference HIRLAM gives the best results. In all experiments the rms error is larger (5-6 ms −1 ) over Iceland than over Scandinavian or European area (3-4 ms −1 ). The behaviour of the 11 km and the 33 km experiments is similar, with the exception of somewhat better pmsl scores in the fine resolution. Figure 8 shows the difference between observed and predicted pmsl and V 10m given by the +48h forecasts of the experiments R33, O33 and B33 over the Icelandic synoptic stations (ca. 10 stations available for the operational data assimilation each night 3 ) classified according to the observed wind velocity. It can be seen that over Iceland, the negative wind bias is due to underprediction of the moderate and strong winds. The error is partly compensated by overprediction of weak winds, especially in the case of the reference experiment R33. The increased drag of B33 leads to retardation of all winds and consequently, to worse overall verification scores. In O33, retardation of the strong winds is the smallest.
Thus, when weak winds dominate (over the whole Scandinavian/European integration area, not shown), the results show a positive bias. In Iceland, moderate and strong winds dominate, leading to a negative wind bias and a positive pressure bias. To correct the systematic error of the low-level wind prediction, weak winds should be retarded and strong winds accelerated in both cases. The proposed orography-related parametrizations do not solve this problem, although the experiments 'O' seem to take a step in the right direction. This is shown in Fig. 9, where the distribution of the lowest model level wind (level 40 centred at the height of 32 m) over Iceland is compared between the experiments O33, B33 and R33. It is seen that the weak winds are weaker and strong winds stronger in O33 compared to R33, while all winds are weaker in B33 compared to R33 and O33. However, the results of R33 and O33 remain quite close to each other, indicating that the improvement obtained by modifying the orography-related surface drag is still insufficient.

Conclusions
Three different subgrid-scale processes contribute to the momentum budget of HIRLAM in the planetary boundary layer (PBL): turbulent diffusion, orographic turbulence related to the SSO and wave-related effects (in fact, blocking of the low-level flow by mesoscale mountains) due to MSO. In the reference model, orographic roughness (z o,oro ) is used to parametrize all orography-related processes. In this study, z o,oro was replaced by new parametrizations for SSO and MSO effects. For the SSO parametrizations, two different approaches were tried: the HIRLAM formulation (eq. 5) and the ECMWF formulation (eq. 6). Needed orographic parameters were derived from fine-resolution ( x = 1 km) digital elevation data.
In order to compare the parametrizations, numerical experiments over Europe were run in January 2000. Coarse (with a horizontal gridlength of 33 km) and fine-resolution (11 km) results were compared. A simple verification of forecasts against (SYNOP) observations turned out to be too smoothing to reveal and understand the differences between the experiments. Momentum fluxes, kinetic energy budget and detailed analysis of observation-forecast differences were thus used. Particular attention was paid to the results over the complex orography of Iceland. Analysis of the surface momentum fluxes showed, as expected, that the influence of the MSO parametrizations decreases when the resolution improves. The surface drag due to the ECMWF-style SSO parametrizations was found to be significantly larger than that given by the reference HIRLAM or the HIRLAM-style SSO parametrizations. It also influenced in deeper layer, dominating over the mesoscale effects even in the coarse resolution experiments. Analysis of the kinetic energy budgets showed that the different parametrizations tend to com-pensate each other. Although the MSO and SSO parametrizations feel the details of different subgrid scales of the underlying orography, the turbulence parametrizations and model dynamics are only able to separate processes larger or smaller than the grid scale, thus possibly smoothing all gradients irrespective of to their origin.
Analysis of forecast-observation differences over different areas showed that the orography-related parametrizations influence in quite local scale. Over flat areas, only minor differences between the experiments were seen. Both the reference HIRLAM and the modified versions tend to overestimate weak and underestimate strong winds in coarse and fine-resolution experiments. A model version with the new HIRLAM-style SSO parametrizations combined with new MSO parametrizations was able to correct this distribution slightly, while the use of the ECMWF-style SSO parametrizations led to a retardation of all winds.
In the future, more extensive model comparisons would be needed in order to fine-tune the new parametrizations. Dependencies between the orography-related parametrizations and the TKE-based turbulence scheme might be analysed further in order to improve the overall result. The behaviour and need of the parametrizations for subgrid-scale orography effects in the fine-scale models (with a kilometre-scale horizontal gridlength and non-hydrostatic dynamics) requires further studies. In this case, even the 1-km resolution elevation data would be insufficiently detailed for derivation of the needed subgrid-scale orographic variables. The diagnostic tools developed and applied in the present study are believed to be useful also in the analysis of the future experiments.

Acknowledgments
Thanks to Carl Fortelius for providing tools for two-dimensional spectral smoothing and for useful comments to the manuscript. The critical comments of Jean-Francois Geleyn and two anonymous reviewers helped to improve the contents and formulations of the article.