The role of horizontal model grid resolution on the development of the daytime boundary layer over mountainous terrain is studied. A simple idealized valley topography with a cross-valley width of 20 km, a valley depth of 1.5 km, and a constant surface heat flux forcing is used to generate upslope flows in a warming valley boundary layer. The goal of this study is to investigate differences in the boundary layer structure of the valley when its topography is either fully resolved, smoothed, or not resolved by the numerical model. This is done by performing both large-eddy (LES) and kilometer-scale simulations with horizontal mesh sizes of 50, 1000, 2000, 4000, 5000, and 10 000 m. In LES mode a valley inversion layer develops, which separates two vertically stacked circulation cells in an upper and lower boundary layer. These structures weaken with decreasing horizontal model grid resolution and change to a convective boundary layer over an elevated plain when the valley is no longer resolved. Mean profiles of the LES run, which are obtained by horizontal averaging over the valley show a three-layer thermal structure and a secondary heat flux maximum at ridge height. Strong smoothing of the valley topography prevents the development of a valley inversion layer with stacked circulation cells and leads to higher valley temperatures due to smaller valley volumes. Additional LES and “1 km” runs over corresponding smoothed valleys reveal that differences occur mainly because of unresolved topography and not because of unresolved turbulence processes. Furthermore, the deactivation of horizontal diffusion improved simulations with 1- and 2-km horizontal resolution.
Under weak synoptic conditions the development of the atmospheric boundary layer is primarily determined by interactions with the earth’s surface. Beside different soil and land-use properties especially the geometry of the land topography influences near-surface atmospheric flows. It is well known that convective boundary layers, which develop over flat plains differ significantly from boundary layers, which evolve over complex terrain (Whiteman 2000). The existence of slopes, valleys, and mountains leads to the generation of thermally driven slope flows, along-valley winds, and plain-to-mountain circulations (Whiteman 2000), which do not exist over homogeneous terrain. Because of their diurnal changing flow patterns and their importance for regional weather in mountainous regions these mesoscale flow regimes have already been intensively studied in the past by means of analytical investigations (e.g., Prandtl 1952; Vergeiner and Dreiseitl 1987), field campaigns (e.g., Weigel and Rotach 2004; Rotach et al. 2008), and modeling studies (e.g., Schumann 1990; Schmidli et al. 2011; Weigel et al. 2007).
Thermally driven flows over mountainous terrain are, however, not only a key feature of the valley boundary layer but also contribute intensely to vertical transport of properties like heat, moisture, momentum, mass, and pollutants (e.g., Gohm et al. 2009). Weigel et al. (2007) showed that moisture and mass fluxes can be much larger over mountainous terrain compared to pure turbulent fluxes over a flat plain. This also affects the state of the free atmosphere, as momentum, heat, and moisture can be transported more effectively to higher altitudes over complex terrain than over flat plains. This was shown by Henne et al. (2004) who measured exchange processes over complex terrain and identified the major role of thermally driven and especially upslope flows in transporting pollutants from the ground to the free atmosphere.
As already stated in Rotach and Zardi (2007) the simulation of thermally driven flows is a real problem for operational numerical weather prediction and climate models, whose horizontal mesh size is still too coarse to properly resolve extremely mountainous topographies such as narrow Alpine valleys. This means that these models either do not capture valleys and mountain slopes with associated valley and slope winds at all or strongly smooth the valley topographies, which can then lead to very different flow patterns compared to reality. An additional source of errors is introduced as coarse-resolution models use boundary layer, surface layer, and turbulence parameterizations, which have been generally developed for conditions over flat terrain (Weigel et al. 2007).
To improve weather and climate models over mountainous regions it is necessary to adapt these parameterization schemes by including effects of unresolved thermally driven flows. In a first step toward this direction the model error in atmospheric boundary layer properties has to be estimated, which results from insufficient resolution of terrain properties and the usage of boundary layer parameterizations developed for flat terrain. The answer to the question if either the unresolved topography or the unresolved turbulence processes are responsible for errors in the boundary layer structure could be very helpful for the development of possible parameterization schemes.
First investigations of slope winds and related heat flux profiles were performed by Noppel and Fiedler (2002). They developed a conceptual model to determine horizontally averaged heat flux profiles in a valley by means of temperature deviation at the slopes, the slope wind layer depth, the stratification in the valley, and the geometrical shape of the valley. It was found that heat transport by slope winds can be larger than turbulent heat fluxes and has a significant impact on the valley boundary layer. Therefore, a representation of slope and valley wind effects in coarse numerical models, which do not resolve valleys, is necessary.
Existing modeling studies of valley and slope winds usually used models with horizontal grid sizes between about 50 m (e.g., Catalano and Moeng 2010; Schmidli 2013) and a few kilometers (e.g., Rampanelli et al. 2004) to study mechanisms that drive slope and valley winds. While slope winds occur due to buoyancy forces associated with horizontal temperature differences between heated or cooled air over slopes and adjacent unperturbed air in the valley center, along-valley winds are the result of horizontal temperature differences between valley and plain boundary layers. The reason for the stronger heating of the valley atmosphere compared to a plain can be explained by different mechanisms: adiabatic subsidence heating, diabatic convective heating, and the valley volume effect. The former suggests that the slope wind circulation induces ascending air masses along the valley slopes and due to mass continuity descending air above the valley center. As the atmosphere above the valley is generally stably stratified, subsidence implies heating of the valley atmosphere (Rampanelli et al. 2004). In addition, a convective boundary layer develops at the valley floor, which contributes to heating of the valley atmosphere from below (Serafin and Zardi 2010). The concept of the valley volume effect is based on the smaller air volume of a valley compared to a plain, which results in stronger heating of the valley atmosphere for an equivalent energy input (e.g., by radiation; Wagner 1938).
In his recent study, Schmidli (2013) shows by means of large-eddy simulations (LES) that subsidence heating is only a local view of perspective and cannot explain the heating of the whole valley atmosphere. According to Schmidli (2013) the energy for the warming of the valley center due to subsidence heating comes from the heated slopes and not from the free atmosphere. Therefore, cross-valley circulations only redistribute heat within the valley and cannot explain the stronger heating of a valley compared to a plain. An appropriate way to explain warmer temperatures in a valley is the valley volume effect, which marks the upper heating rate limit of the valley boundary layer (Schmidli 2013).
In this study we are interested in systematically investigating the impact of grid resolution on the resulting structure of thermally induced flows and the overall exchange between the valley atmosphere and the free troposphere aloft. To our knowledge similar investigations are missing in the current research. For this purpose we gradually decrease the grid size from high-resolution LES to coarse resolution at the lower end of hydrostatic modeling. The differences in the flow and turbulence structure are assessed between simulations with a fully resolved valley and with various degrees of smoothing of the topography including a completely flat topography. Additional runs in LES mode and with 1-km horizontal mesh size over corresponding smoothed topographies are performed to distinguish between effects of valley smoothing and unresolved turbulence structures. Furthermore, the impact of horizontal diffusion on the flow structure in kilometer-scale simulations is investigated by performing sensitivity runs with deactivated horizontal smoothing. To reduce complexity of the developing flow, a very simple straight valley without an adjacent plain is chosen in this investigation. A similar setup was applied in a recent study of Schmidli (2013). Following the methods of this paper, the comparison of LES and kilometer-scale simulations in this study is enabled by temporal and spatial averaging of the flow fields. This also allows us to split fluxes into mean advective, resolved turbulent, and subgrid-scale components, which helps to clarify differences between coarse- and high-resolution model runs.
The paper is organized as follows: in section 2 the numerical model and the experimental setup is described. Section 3 explains the averaging method, the flow decomposition into mean advective, resolved turbulent and subgrid-scale parts, and the computation of mean vertical profiles. Section 4 exhibits the simulation results and differences in boundary layer structures due to varying horizontal grid sizes. The conclusions are given in section 5.
2. Model and setup
In this study the Advanced Research version of the Weather Research and Forecasting Model (ARW-WRF), version 3.4 is used for numerical simulations (Skamarock et al. 2008). WRF has already been succesfully applied for idealized simulations of thermally driven flows in the kilometer scale (Rampanelli et al. 2004; Schmidli et al. 2011) and for LES studies (Catalano and Moeng 2010; Catalano and Cenedese 2010).
WRF is a nonhydrostatic, fully compressible numerical model, which uses a horizontally staggered Arakawa C grid with a terrain-following vertical coordinate based on the dry-hydrostatic pressure (Skamarock et al. 2008). In this study a third-order Runge–Kutta time integration scheme, fifth-order horizontal, and third-order vertical advection scheme is used. In LES mode, subgrid-scale turbulence is parameterized by a three-dimensional 1.5-order turbulent kinetic energy (TKE) closure (Deardorff 1980). For kilometer-scale simulations subgrid-scale vertical mixing is done by a Mellor–Yamada–Nakanishi–Niino (MYNN) level-2.5 boundary layer scheme (Nakanishi and Niino 2006) and a Smagorinsky first-order closure with a Smagorinsky constant of Cs = 0.25 (Skamarock et al. 2008) is used for horizontal diffusion in physical space (see Skamarock et al. 2008). The MYNN scheme is used, as it is a local TKE scheme, which was already successfully used in WRF simulations of Schmidli et al. (2011) and provides better results over complex terrain than the Yonsei University (YSU) boundary layer parameterization scheme (Hong et al. 2006; Schmidli et al. 2011, see section 4). At the surface a Monin–Obukhov similarity scheme (Monin and Obukhov 1954) using four stability regimes of Zhang and Anthes (1982) is applied. The decomposition of the turbulent flow into resolved and mean components is done according to the method described in section 3. To reduce the amount of data storage needed for computations, a statistics module has been implemented in WRF, which allows for an online averaging and flux computation while the model is integrating.
The model terrain is based on a simple infinitely long and straight valley topography, which is identical to the one used in Schmidli (2013). The modeling domain has an extension of 40 km in the cross valley and 10 km in the along-valley direction. The straight valley is defined by two parallel mountain ridges with 1.5-km height and a crest-to-crest distance L of 20 km. Both the valley floor and the mountain ridges consist of flat areas with a width of 1 km in cross-valley direction (see Fig. 1). The maximum inclination of the mountain slopes is 14.7°. As no adjacent plain exists in the topography, only slope winds and no along-valley winds develop in this setup.
LES runs use horizontal mesh sizes of 50 m and 175 vertically stretched levels with varying depths of 8 m close to the ground increasing to 50 m higher aloft. For kilometer-scale runs horizontal grid distances of 1, 2, 4, 5, and 10 km and a vertically stretched grid with 60 levels and level distances of 25–100 m are used. In the LES setup the lowest full model level is located at 4 m, whereas in the kilometer-scale runs it is positioned at about 12.5 m above ground level (AGL). The time step is 0.5 and 6 s for LES and kilometer-scale runs, respectively. Test runs with larger time steps for kilometer-scale runs do not change the results compared to a 6-s time step.
The model top is set to 5 km with a Rayleigh damping layer covering the uppermost 1000 m. Sensitivity runs with higher model tops at 8 and 20 km and correspondingly deeper damping layers of 3 and 10 km, respectively, did not show any significant gravity waves above 3-km height and, hence, no differences in the boundary layer structure compared to the runs with a 5-km model top. In the horizontal direction periodic lateral boundary conditions are applied in both cross- and along-valley directions to simulate infinitely long and repeating parallel valleys.
The simulations are initialized with an atmosphere at rest, a potential temperature of 297 K at sea level height, and a constant potential temperature gradient of 3 K km−1. To avoid moist processes, a constant relative humidity of 40% is chosen at the beginning of the simulations. The thermal forcing is defined by a constant surface sensible heat flux of 150 W m−2 and the surface roughness length is set to 0.16 m. To trigger convection at the beginning of the simulations, randomly distributed potential temperature perturbations with a maximum amplitude of 0.5 K are added to the five lowermost model levels. All simulations are run for 5 h with a data output interval of 1 h.
The adaptation of the valley shape to the varying horizontal model resolutions is done by interpolating a highly resolved topography to the respective coarse model grid and smoothing the resulting terrain with the aid of the smooth–desmooth operator of Shapiro (1970):
Here Zi,j is the interpolated two-dimensional topographic height field at coarse model grid points i and j, is the corresponding smoothed field, and S is a smoothing coefficient. A single pass of the operator consists of a smoothing step (S = 0.50) to remove two-grid-interval waves (Shapiro 1970) and smooth short wavelengths. It is followed by a desmoothing step (S = −0.52) to restore original waves (Grell et al. 1994). This smoothing method is adopted from the WRF real-case preprocessing system (WPS; Wang et al. 2012), where it is applied to smooth high-resolution terrestrial data.
The smoothing leads to a reduction of the mountain heights and to elevated valley floors for coarser model runs. The relative resolution Δx/L with horizontal mesh size Δx and valley width L = 20 km can be used to define the fully resolved (Δx/L ≤ 0.1) and smoothed (Δx/L > 0.1) valley topography (see Table 1). In the extreme case of a simulation with 10-km horizontal mesh size, the valley topography becomes an elevated plain at 750-m height. The two reference cases used in this study are LES runs with 50-m horizontal mesh size, where the first one uses the fully resolved valley topography (dx50m) and the second one applies an elevated plain (dx50mT10) to assess the simulation with 10-km grid spacing. The other kilometer-scale runs with 1- and 2-km grid spacing resolve the valley well (Δx/L=0.05 and 0.1, respectively), whereas the cases with 4- and 5-km horizontal resolution result in strong valley smoothing (Δx/L = 0.2 and 0.25, respectively). To investigate the role of unresolved turbulence motions, additional runs in the LES mode and with a 1-km horizontal mesh size are performed over corresponding smoothed topographies. The type of topography is indicated by respective suffixes in the simulation names (e.g., “T5” for dx5km topography). Additional kilometer-scale runs with deactivated horizontal diffusion schemes are conducted to study the role of horizontal diffusion on the flow structure. The deactivation of horizontal diffusion is indicated by the suffix “ND” in simulation names. An overview of the different simulations and their characteristics is given in Table 1.
3. Averaging method
a. Flow decomposition
To enable the comparison of LES and kilometer-scale simulations the model output has to be decomposed into mean advective, resolved turbulent, and subgrid-scale parts. Therefore, the flow variables are averaged following the approach of Schmidli (2013). The fully turbulent variable is split into a model gridbox average , which is dependent on mesh size and time step of the numerical model and a subgrid-scale part ψ′(x, t):
The resolved turbulent part ψ″(x, t) is then computed from the gridbox average by
where the time and along-valley averaging operator 〈 〉 is defined as
with the time-averaging interval T = 40 min and the along-valley averaging interval Ly = 10 km. The chosen averaging intervals are applied to be in agreement with the settings of Schmidli (2013), whose simulations are used to test the reliability of our results (see section 4a). Time averaging is based on a sample interval of 1 min.
b. Computation of vertical fluxes
By averaging over a model grid box the purely turbulent vertical flux becomes and can be split into a gridbox average (GBA) and a subgrid-scale (SGS) part by means of Eq. (2) as
where it is assumed that and . The GBA part can be further decomposed into a mean advective (MEAN) and resolved turbulent part (RES) by the usage of Eq. (3) and the averaging operator defined in Eq. (4):
with the assumption that and .
The MEAN and RES fluxes can be calculated directly from the model output, whereas the subgrid-scale fluxes are computed by means of diffusivity coefficients and corresponding gradients of gridbox-averaged variables:
where KHV is a vertical diffusivity coefficient for scale quantities, such as heat. The diffusivity coefficients are computed by the respective turbulence model (3D TKE Deardorff scheme for LES and MYNN level-2.5 model for kilometer-scale runs).
A special procedure has to be performed for the computation of the mean vertical heat flux . According to thermodynamic theory (Bohren and Albrecht 1998) only the change in enthalpy dh = cpdθ, with enthalpy h, specific heat capacity for constant pressure cp, and potential temperature θ is defined exactly, but not its integral h = cpθ, which is dependent on the definition of a reference state. Therefore, only fluxes in the form are unique. In this study is defined according to Schmidli (2013) as , where θ0(z) is the spatial average of at a certain height level.
To get mean vertical profiles over the whole valley region, the fluxes in Eq. (6), as well as potential temperature and heating rate are averaged along constant height levels within and above the valley volume (x = −10 to x = 10 km). This operation requires a data interpolation from model levels to Cartesian coordinates. The averaging levels have a vertical distance of 20 m and are illustrated schematically with horizontal lines in Fig. 1.
4. Simulation results
a. Instantaneous flow
In this section an overview of the instantaneous flow field after 5 h of simulation is given for runs with different horizontal model resolutions. As the model setup is similar to the Advanced Regional Prediction System (ARPS) simulations of Schmidli (2013), WRF-LES results (dx50m) were compared to these data. This comparison reveals a very good agreement between both models (not shown) and enables to use the dx50m simulation as a reference case.
As in Schmidli (2013) the atmospheric flow (dx50m) is fully turbulent after about 2 h of simulation. In Fig. 2 snapshots of cross-valley and vertical wind speed after 5 h of integration are shown. Upslope winds are well defined with upslope wind speeds between 2.2 and 3.5 m s−1 in the valley-resolving simulations (dx50m and dx1km), but are weaker (up to 1.3 m s−1) in the very shallow valley of the dx5km case (Fig. 2). In the dx50m run two cross-valley circulation cells on top of each other with upward motion near the slopes and horizontal return flows toward the valley center above and below the mountain ridges can be identified. The stacked cells are weakly developed in the dx1km simulation especially in the lower valley atmosphere. In the dx4km case, only one big circulation cell exists with a much deeper upslope-wind layer and slope-parallel return flows. Because of the lack of a valley topography and the coarse-resolution negligible horizontal winds evolve in the convective boundary layer over the elevated plain in the dx10km case.
Horizontal cross sections of vertical wind speed at about 100 m AGL show convective thermals, which are advected upslope in the dx50m case (Fig. 2). These cells show typical hexagonal structures with horizontal diameters of about 1–2.5 km. These cells are very similar to thermals evolving in a convective boundary layer over flat terrain (e.g., Schmidt and Schumann 1989) and occur also in the LES of the elevated plain (dx50mT10, not shown). With increasing horizontal mesh sizes these finescale thermals lose their structure (dx1km) and finally vanish in the dx4km, dx5km, and dx10km runs.
b. Averaged flow
The spatial and temporal averaging of the flow according to the method described in section 3 leads to more comparable fields. In Fig. 3, cross sections after 5 h of simulation are shown for runs with different horizontal resolutions. In Fig. 4, runs with mesh sizes of 1 km are displayed over differently smoothed valleys (e.g., dx1kmT4), as well as the 2-km mesh size run with deactivated horizontal diffusion (dx2kmND). Two vertically stacked circulation cells are clearly visible in the dx50m case, but are only very weakly developed in the simulation dx1km (Fig. 3). The cells are separated by a valley inversion layer, which is located below mountain ridge height. The valley inversion is still weakly present in the dx2km case, but missing in the dx4km and dx5km runs. In line with the disappearance of this valley inversion only one big cross-valley circulation develops with slope-parallel return flows and strong subsidence in the valley center (e.g., dx2km, dx4km, and dx5km; Fig. 3). In the dx50m and dx1km runs the return flows are horizontally oriented and subsidence over the valley center is divided into two regions of maximum descent with values as low as −0.47 m s−1 above and below the valley inversion. It can also be seen that in the dx4km case the upslope wind layer becomes deeper and maximum upslope wind speeds (up to 2.1 m s−1) occur at an altitude of about half the valley depth, whereas maximum slope winds (up to 2.2 m s−1) are found near the mountain tops in the dx50m and dx1km runs. In contrast, upslope winds and subsidence over the valley center are comparatively weak with values of 1.3 and −0.25 m s−1, respectively, in the shallow valley of the dx5km case.
The convergence zones over the mountain ridges are characterized by convective cells with mean updrafts of up to 1.8 m s−1 in the dx50m case. These narrow thermals become weaker and more expanded in the kilometer-scale runs with maximum updrafts of only 0.4 and 0.3 m s−1 in the dx4km and dx5km case, respectively. The larger width of the cells in the dx4km case is primarily caused by coarser horizontal resolution.
The comparison of the cross sections in Fig. 3 makes it clear that the boundary layer structure cannot be described as one mixed layer in the dx50m and dx1km cases. In fact, the boundary layer has to be divided into a lower valley boundary layer and an upper weakly stable layer, which are separated by the valley inversion layer. Because of these two layers there are also two boundary layer heights: a lower one, which marks the height of the first inversion and an upper one, which marks the transition to the free atmosphere. Observations of pollution concentrations in valleys show that the lower valley inversion also exists in reality and often traps certain properties in the lower part of the valley atmosphere (Henne et al. 2004; Gohm et al. 2009).
Because of the existence of two mixed layers in the dx50m and dx1km runs, two boundary layer heights are introduced according to the definition of Catalano and Moeng (2010). The lower boundary layer height is determined as the altitude where the potential temperature gradient exceeds the value of 0.001 K m−1 when moving upward from the surface. Correspondingly, the upper boundary layer height is located at the altitude where the gradient falls below this threshold when moving downward from the model top.
In the dx50m and dx1km runs a clear separation of the upper and lower boundary layer height is visible (Fig. 3). The difference of the two heights is largest over the valley center, with a lower boundary layer height below the mountain crests. Over the ridges both heights are almost equal. In the dx4km and dx5km runs there is no difference between the upper and lower boundary layer height due to the lack of a valley inversion.
To separate the effect of resolved topography and resolved turbulent motions all smoothed topographies are used for runs in LES mode and with mesh sizes of 1 km (Fig. 4). Compared to the dx2km run the dx1kmT2 case is able to simulate the vertically stacked circulation cells and a weak valley inversion layer. The additional dx2kmND run with a deactivated horizontal diffusion scheme is able to resolve both the circulation cells and the valley inversion and agrees much better with the reference case than the dx2km run (cf. Fig. 3). This means that horizontal smoothing by the diffusion scheme might be too strong in kilometer-scale simulations. Comparisons of dx1kmT4 and dx1kmT5 with dx4km and dx5km cases, respectively, show that in all runs only one circulation cell and no valley inversion exists. Because of the finer horizontal resolution the upslope wind maxima are shifted toward the mountain ridges and updrafts over the crests are narrower and stronger in the dx1kmT4 and dx1kmT5 cases compared to the coarser-resolution runs (Figs. 3 and 4). The dx1kmT10 run is able to resolve smaller convective cells, but is very similar to the dx10km run (not shown) in terms of thermal structure of the boundary layer.
Mean vertical profiles of cross-valley wind speed after 5 h of simulation are shown in Fig. 5 for different cross-valley locations. Note that a linear interpolation method is used to derive profiles at equal cross-valley positions for all simulations. To better distinguish the large number of simulations from each other, the following convention has been introduced and applied to all line graphs: the color indicates the type of topography, whereas the line style specifies the simulation type (see Table 1).
Upslope flows near the surface increase from about 0.9 m s−1 at the position x = 2 km to about 2.1 m s−1 near the mountain ridge at x = 8 km in the dx50m simulation (Fig. 5). The agreement between kilometer-scale runs and the LES (dx50m) is best near the mountain top at x = 8 km (Fig. 5c) as the near-surface flow convergence and the return flow toward the valley center is similar in all valley resolving simulations near the mountain ridges (see Figs. 3 and 4). However, in the valley center and at midslope positions (x = 2 and x = 4 km) the discrepancy between different simulations in the upslope wind profile is large due to unresolved vertically stacked circulation cells, especially in the dx4km and dx5km cases. The stacked circulation cells are clearly visible in the dx50m simulation at x = 2 and x = 4 km by means of the alternating positive and negative values of the cross-valley wind component with height. This vertical flow structure is also captured by the dx1km run. The deactivation of horizontal diffusion in the dx2kmND run leads to much better results than in the dx2km case, as vertically stacked circulation cells are now visible in the vertical profile at x = 2 km (Fig. 5a) and as the profiles are closer to the reference run (dx50m). Differences of upslope winds between dx1kmT4, dx1kmT5 and dx4km, dx5km cases, respectively, are largest near the mountain ridges, as coarse-resolution runs simulate strongest upslope winds at midslope positions (cf. Fig. 4).
Upslope winds near the surface show a sharp jetlike structure in the dx50m run with a wind speed maximum at about 20 m (x = 2 and x = 4 km) and 36 m AGL (x = 8 km). This jet structure is not captured by kilometer-scale simulations. An additional run of the dx50m case, which uses the coarser vertical model grid of the kilometer-scale runs does not properly resolve the shallow surface jet either (not shown). This indicates the importance of a high vertical grid resolution near the surface. In the elevated plain simulations (dx10km, dx1kmT10, and dx50mT10) a convective boundary layer develops without any mean directed flow.
c. Spatially averaged profiles
To study differences in the boundary layer structure, which result when the valley is fully resolved, smoothed, or unresolved, horizontally averaged profiles over the valley region are calculated for potential temperature, vertical sensible heat flux, and heating rate.
The time evolution of mean potential temperature profiles is displayed in Fig. 6. The dx1km run agrees very well with the reference case and is able to simulate the valley inversion layer. The dx2km run simulates too warm temperatures in the lower boundary layer and does not capture the inversion layer. The deactivation of horizontal diffusion for 1- and 2-km mesh sizes (dx1kmND, dx2kmND) significantly improves the results and enables the formation of the valley inversion layer for dx2kmND and a nearly perfect agreement with the reference case for dx1kmND. When the valley is smoothed (Δx/L ≥ 0.2) temperatures are much too warm, especially at the beginning of the simulations. In general, the development of the thermal structure is very similar for runs with identical topographies but different horizontal mesh sizes (e.g., dx4km, dx1kmT4, dx50mT4), whereas large differences occur in runs with identical mesh sizes but with different topographies (e.g., dx50m, dx50mT4, dx50mT5). This implies that the largest differences in the thermal structure of the boundary layer occur due to unresolved or smoothed valley topographies and not due to unresolved turbulence motions. In addition, the time evolution of potential temperature profiles is very similar for runs with identical topographies. For example the growth rate of the boundary layer is similar in the dx50mT10, dx1kmT10, and dx10km cases as can be seen in the time development of potential temperature profiles (Fig. 6).
Potential temperature profiles of the dx50m and dx1km runs develop a so-called three-layer thermal structure (e.g., Vergeiner and Dreiseitl 1987; Schmidli 2013) with a well-mixed lower valley boundary layer, a valley inversion layer, and an upper weakly stable layer. The valley inversion layer develops due to a local heating rate minimum. This means that stronger heating above the region of the developing inversion leads to an increasing potential temperature gradient in between.
This local heating rate minimum is clearly visible in vertical heating rate profiles over the valley center (at x = 0 km) after 2 h of simulation in the dx50m and dx1km runs (Fig. 7). The total heating rate (TOT) is split according to Eq. (4) in Schmidli (2013) into a heating rate due to mean advection (ADV) and into resolved and subgrid-scale turbulent heat flux divergence (TRB). Near the surface the atmosphere is heated as a result of turbulent heat flux divergence [bottom-up heating as discussed in Serafin and Zardi (2010)]. This leads to the development of a convective boundary layer with rising thermals. The turbulent heating rate (TRB) decreases with height and becomes negative (i.e., a cooling rate) in the entrainment layer of the lower boundary layer. This is the region where rising thermals reach the stably stratified part of the valley atmosphere and produce a cooling effect due to the turbulent exchange of the overshooting potentially cooler thermals (cf. Schmidli 2013). This mechanism is visible in both LES and kilometer-scale runs.
Major differences occur in profiles of mean advective subsidence heating (Fig. 7). When the valley is well resolved (dx50m and dx1km runs) strong subsidence heating (ADV) develops below and above the stably stratified part of the valley atmosphere with a local heating rate minimum in between. With decreasing horizontal model resolution this local minimum in subsidence heating rate weakens (dx2km) and finally vanishes (dx4km). In the dx4km case only one peak of subsidence heating remains due to the strongly reduced valley depth.
Horizontal averaging over the valley region (Fig. 8) smooths the profiles compared to the ones at the valley center shown in Fig. 7. However, the local heating rate minimum is still visible in the averaged profiles of the dx50m case between 2 and 4 h (Fig. 8). In the dx1km run the local minimum is no longer visible, but a heating rate decrease occurs at the altitude of the inversion layer. While heating rate profiles are very different during the first 2 and 3 h of simulation between cases with different valley topographies (Figs. 8a,b), they become very similar during the last two hours (Figs. 8c,d) due to the growing convective boundary layer. After 5 h the boundary layer is homogeneously heated in runs over smoothed and unresolved valleys (e.g., dx4km, dx5km, dx10km), whereas strong heating occurs below and weaker heating above the valley inversion layer in runs with fully resolved topographies (e.g., dx50m, dx1km, dx2km).
The different heating rate profiles between the set of simulations with well resolved (dx50m, dx1km, dx2km) and smoothed (dx4km, dx5km) valleys might also result from different valley volumes and different area–height distributions (Steinacker 1984). For example, the valley air volume in the dx4km and dx5km runs is only one-half and one-third of the fully resolved valley volume in the dx50m case, respectively. These much smaller valley volumes consequently support stronger heating, especially at the beginning of the simulations (e.g., Fig. 8a).
To study differences in heating of the boundary layer in dependence of different horizontal grid resolutions, potential temperature tendencies have been density weighted and integrated over the corresponding boundary layer air volume according to Eq. (5) in Schmidli (2013). The air volume is defined from the surface to the upper boundary layer height (cf. Fig. 3), which represents the overall depth affected by heating (see Table 1 for the different valley floor heights). A time series of the resulting heating rates and boundary layer volumes is shown in Fig. 9. When studying runs with different grid sizes, heated volumes are smaller in the dx4km, dx5km, and dx10km runs compared to the dx50m, dx1km, and dx2km runs at the beginning of the simulations due to lower boundary layer heights and elevated terrain surfaces. Because of the smaller air volumes over smoothed topography the boundary layer is stronger heated and potential temperatures at respective heights are higher at the beginning of the simulations (cf. Fig. 6). As time proceeds the heating rate strongly decreases in the dx4km, dx5km, and dx10km runs (Fig. 9) due to increasing boundary layer heights and corresponding increasing boundary layer air volumes (Fig. 9). Runs with 1-km mesh size over smoothed topography (e.g., dx1kmT2, dx1kmT4, dx1kmT5) reveal larger volumes than corresponding coarse-resolution runs (dx2km, dx4km, dx5km). This difference is due to deeper convective cells over the mountain ridges in the 1-km runs (cf. Figs. 3 and 4). In the dx50m, dx1km, and dx2km runs, the upper boundary layer height and consequently the heated air volume grow slower with time, as the upper boundary layer height is located much higher at the beginning of the simulations due to higher mountains (cf. Fig. 6). This results in a weaker decrease of volume-averaged heating rate (Fig. 9). After 5-h integration time the heating rates of simulations with strongly smoothed topographies approach the dx50m, dx1km, and dx2km runs, but are still larger because of smaller volumes. Surprisingly, the heating rates of the dx10km and dx50mT10 simulations over the elevated plain are very similar, although the boundary layer volume of the dx50mT10 run is considerably larger. The reason for similar heating rates despite of different volumes is a deeper boundary layer with stronger entrainment near the boundary layer top (not shown) in the dx50mT10 run compared to the dx10km run (cf. Fig. 6). The additional upper part of the volume in the dx50mT10 run is heated by downward mixing of potentially warmer air from above (cf. Fig. 8d). Hence, the MYNN boundary layer scheme is not able to correctly represent the entrainment process and consequently the boundary layer depth. Additional runs of the kilometer-scale cases with the YSU boundary layer parameterization scheme (Hong et al. 2006) agree perfectly with the LES run over the elevated plain (dx10km) in terms of boundary layer air volume, vertical profiles of potential temperature, and heating rate (not shown). Over the valley terrain, however, the YSU dx1km, dx2km, dx4km, and dx5km runs do not capture the upslope winds and thermal structure properly and tend to develop deeper boundary layers than in the runs with the MYNN scheme. Similar results were found in the model comparison study of Schmidli et al. (2011). This shows that the strength and structure of thermally driven flows are very sensitive to the choice of the boundary layer parameterization.
Figure 10 shows valley mean profiles of vertical sensible heat flux, which is split into mean advective, resolved turbulent, and subgrid-scale parts according to the method described in section 3b. Because of the development of upslope circulations, all valley resolving simulations (e.g., dx50m, dx1km, dx2km) show nonzero mean advective heat flux profiles with maximum fluxes at heights with strongest upslope winds. Over the elevated plains, however, mean advective vertical heat fluxes are very small. Because of the broad convective cells over the wide mountain plateaus the dx4km case (cf. Fig. 3) shows the greatest deviation in mean vertical heat flux from the reference run (dx50m).
Resolved turbulent heat fluxes are strongest in the convective boundary layers over the valley floor (Fig. 10b). LES runs over well-resolved valleys show secondary resolved turbulent heat flux maxima over the mountain ridges (e.g., dx50m, dx50mT1, dx50mT2). Over smoothed valleys (dx50mT4, dx50mT5) this secondary maxima becomes weaker and finally vanishes over the elevated plateau (dx50mT10). Subgrid-scale heat fluxes are more dominant in kilometer-scale runs with largest values at the valley floor (Fig. 10c). While over the elevated plains the subgrid-scale heat flux decreases nearly linearly with height, the valley resolving simulations indicate a local heat flux maximum over the mountain ridges, which act as an additional elevated heating surface for the atmosphere.
The combination of the three parts results in mean profiles of total vertical heat flux (Fig. 10d). Surprisingly the strength of the linear decrease of the heat flux in the valley is similar to the profile over the elevated plain. At mountain crest height a secondary heat flux maximum is visible, with a linear heat flux decrease above. This secondary maximum is a result of the elevated heating area of the ridges and the converging upslope winds, which induce turbulence and strong mean updrafts. Total heat fluxes are very similar for runs with identical topographies (e.g.,dx2km, dx1kmT2, dx50mT2).
d. Evolution of boundary layer height
One main difference between the simulations with variable horizontal grid resolutions lies in the capability of developing a valley inversion layer, which divides the boundary layer into an upper and lower part. The upper and lower heights of these two subboundary layers are calculated as described in section 4b. Figure 11 shows the evolution of valley averaged potential temperature profiles and corresponding upper (dashed lines) and lower (dotted lines) boundary layer heights for dx50m, dx2km, and dx4km runs. Profiles for the dx1km run are not shown, as they are very close to the dx50m case. While in the dx50m run the valley inversion layer still exists after 5 h of simulation, it has disappeared in the dx4km and dx2km runs after approximately 2 and 3 h, respectively. After the decay of the inversion layer in the dx2km and dx4km runs a weakly stable layer remains with an upper boundary layer height at about 2 km. The evolution of the upper boundary layer height is strongly dependent on the heights of the mountains. Over well-resolved topographies the upper boundary layer height grows only slowly, whereas a strong increase of the boundary layer height is visible over strongly smoothed terrain (Fig. 11). However, after 5 h of integration the upper boundary layer heights of the different simulations converge and are only about 300 and 120 m lower in the dx2km and dx4km runs, respectively, compared to the dx50m simulation.
To study the influence of the heat flux forcing on the boundary layer evolution, simulations are performed with different surface sensible heat fluxes of 50, 100, and 200 W m−2. The corresponding valley averaged potential temperature profiles after 5 h of simulation are shown in Fig. 12. In the dx50m and dx2km runs valley inversion layers exist for heat flux forcings up to 150 and 100 W m−2, respectively. For stronger surface heat fluxes the valley inversion is eroded due to bottom-up and subsidence heating and a weakly stable mixed layer with only one circulation cell remains. For strongly smoothed valleys (dx4km) only an upper boundary layer height exists after 5 h independently of the chosen surface heat flux. The difference in the upper boundary layer heights is largest for small heat flux forcings and reduces to about 150 m for a heat flux of 200 W m−2.
In this study, the influence of the horizontal model grid resolution on the boundary layer structure over an idealized valley was investigated. The goal was to study differences of atmospheric flows over complex terrain, which occur when the topography is fully resolved, smoothed, or unresolved in a numerical model. For this purpose, both LES and kilometer-scale simulations with horizontal grid distances of 50 m, as well as 1, 2, 4, 5, and 10 km were performed with WRF. The valley was fully resolved in the dx50m, dx1km, and dx2km runs, smoothed in the dx4km and dx5km runs, and unresolved (i.e., resulting in an elevated plain) in the dx10km run. To separate the influence of resolved topography and resolved turbulent motions additional runs in LES mode and with 1-km horizontal mesh size were performed over corresponding smoothed valleys. In addition, all kilometer-scale cases were rerun with deactivated horizontal diffusion scheme.
Large differences were found in the cross-valley wind field and temperature structure in the valley due to both different horizontal and vertical grid resolutions and different valley geometries as a result of topography smoothing. The correct representation of the valley depth was crucial for an appropriate simulation of the flow. Major differences occurred between a horizontal grid resolution of 2 km (Δx/L = 0.1) and 4 km (Δx/L = 0.2), that is, between a fully resolved and smoothed valley, due to different valley depths. Surprisingly, for the chosen valley, simulations with 1-km horizontal spacing produced very similar flow structures as the large-eddy simulation. Given the horizontal scale of the chosen valley it is concluded that about 10–20 grid points (Δx/L ≤ 0.1) over the entire valley topography are required to appropriately reproduce the flow structure over highly complex terrain. Additionally, it could be demonstrated that horizontal diffusion seemed to be too strong in kilometer-scale runs. Sensitivity runs with horizontal mesh sizes of 1 and 2 km (Δx/L = 0.05 and 0.1, respectively) and deactivated horizontal diffusion scheme revealed results that were much closer to the reference run compared to cases with horizontal diffusion.
The most important characteristics of the simulated flow were the development of a valley inversion and the formation of two vertically stacked circulation cells. Both features were captured by the LES and the dx1km runs. The development of the valley inversion was enabled by a local heating rate minimum with maximum heating above the inversion layer. Similar valley inversions and vertically stacked circulation cells were already found in Serafin and Zardi (2010) and Schmidli (2013). Simulations with smoothed valleys were not able to produce this three-layer thermal structure and showed only one circulation cell.
In addition, it was found that potential temperatures in the boundary layer were too high in simulations with smoothed terrain compared to the reference run. The reason is larger average heating rates in the boundary layers over smoothed valleys due to smaller boundary layer air volumes, especially at the beginning of the simulations. As time proceeded the heated air volumes strongly increased resulting in decreasing heating rates in coarse-resolution runs. In the LES and dx1km cases, boundary layer volumes were larger and increased only slowly with time, which led to lower heating rates and potential temperatures.
Differences between the simulations were also found in vertical profiles of upslope winds. The depth of the upslope wind layers was too large and the near-surface jet structure was not captured in runs with both coarse horizontal resolutions and large vertical level distances near the surface. In addition, subsidence in the valley center was strong when only one circulation cell was present (e.g., dx4km run).
The computation of horizontally averaged profiles of vertical sensible heat flux showed further differences. The typical linear decrease of total heat flux with height over the elevated plain was modified by the valley topography especially at mountain crest height resulting in a secondary heat flux maximum. In addition, values of mean advective vertical heat fluxes were nonzero in the valley due to upslope winds in contrast to situations over elevated plains.
The development of the valley inversion layer in the dx50m, dx1km, dx1kmND, and dx2kmND runs required the definition of a lower and upper boundary layer height. While the lower boundary layer height only existed in runs with Δx/L ≤ 0.1, the upper boundary layer heights were very similar in all simulations after 5 h of integration. It was found that the development of the upper boundary layer height is strongly dependent on the height of the mountains, as upper boundary layer heights strongly increased with time over smoothed topographies, but stayed almost constant over well-resolved terrain with high mountains.
Sensitivity runs with different sensible heat flux forcings showed that differences in the upper boundary layer height were largest for small surface heat fluxes (e.g., 50 W m−2). We suppose that this difference is even larger for deeper and narrower valleys. To test this hypothesis, additional investigations will be necessary, which should also take into account the development of along-valley winds.
The present results show that smoothed or unresolved valley topographies in large-scale models lead to errors in the boundary layer structure, such as different stratification, heating rates, and boundary layer depths (from the surface to the upper boundary layer height) due to deficient or unresolved thermally driven flows. We suppose that these errors can increase with more complex valley topographies, inhomogeneous surface forcings, and the inclusion of moist processes. The additional runs in LES mode and with 1-km spacing over smoothed terrain revealed that major differences in the boundary layer structure occur due to unresolved topography, while unresolved turbulent motions are less important for this valley setup. This means that the ratio Δx/L seems to determine, whether the flow can be simulated properly. In this study an appropriate simulation of the boundary layer structure was able for values Δx/L ≤ 0.1 and a fixed valley width L = 20 km. Further studies would be necessary to generalize this result by varying the valley width. Therefore, future investigations should focus on the impact of different valley geometries (e.g., valley depth, width, curvature), atmospheric background conditions (e.g., stratification, synoptic wind), and land-use types on thermally driven flows and related exchange processes. Results of such studies could serve as a basis for developing boundary layer parameterizations in complex terrain for use in large-scale operational weather prediction and climate models.
This work was supported by the Austrian Science Fund (FWF) under Grant P23918-N21 and by the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Research Platform Scientific Computing at the University of Innsbruck. This study profited from in-depth discussions with Jürg Schmidli from ETHZ and comparisons with his results.