Journal topic
The Cryosphere, 12, 3137–3160, 2018
https://doi.org/10.5194/tc-12-3137-2018
The Cryosphere, 12, 3137–3160, 2018
https://doi.org/10.5194/tc-12-3137-2018

Research article 04 Oct 2018

Research article | 04 Oct 2018

# Spatial variability in snow precipitation and accumulation in COSMO–WRF simulations and radar estimations over complex terrain

Spatial variability in snow precipitation and accumulation in COSMO–WRF simulations and radar estimations over complex terrain
Franziska Gerber1,2, Nikola Besic3,4, Varun Sharma1, Rebecca Mott2,5, Megan Daniels6, Marco Gabella4, Alexis Berne3, Urs Germann4, and Michael Lehning1,2 Franziska Gerber et al.
• 1Laboratory of Cryospheric Sciences, School of Architecture, Civil and Environmental Engineering, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
• 2WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland
• 3Environmental Remote Sensing Laboratory, School of Architecture, Civil and Environmental Engineering, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
• 4Radar, Satellite, Nowcasting Department, MeteoSwiss, Locarno, Switzerland
• 5Institute of Meteorology and Climate Research, Atmospheric Environmental Research (KIT/IMK-IFU), KIT-Campus Alpin, Garmisch-Partenkirchen, Germany
• 6independent researcher: Sydney, Australia

Correspondence: Franziska Gerber (gerberf@slf.ch)

Abstract

Snow distribution in complex alpine terrain and its evolution in the future climate is important in a variety of applications including hydropower, avalanche forecasting and freshwater resources. However, it is still challenging to quantitatively forecast precipitation, especially over complex terrain where the interaction between local wind and precipitation fields strongly affects snow distribution at the mountain ridge scale. Therefore, it is essential to retrieve high-resolution information about precipitation processes over complex terrain. Here, we present very-high-resolution Weather Research and Forecasting model (WRF) simulations (COSMO–WRF), which are initialized by 2.2 km resolution Consortium for Small-scale Modeling (COSMO) analysis. To assess the ability of COSMO–WRF to represent spatial snow precipitation patterns, they are validated against operational weather radar measurements. Estimated COSMO–WRF precipitation is generally higher than estimated radar precipitation, most likely due to an overestimation of orographic precipitation enhancement in the model. The high precipitation amounts also lead to a higher spatial variability in the model compared to radar estimates. Overall, an autocorrelation and scale analysis of radar and COSMO–WRF precipitation patterns at a horizontal grid spacing of 450 m show that COSMO–WRF captures the spatial variability normalized by the domain-wide variability in precipitation patterns down to the scale of a few kilometers. However, simulated precipitation patterns systematically show a lower variability on the smallest scales of a few hundred meters compared to radar estimates. A comparison of spatial variability for different model resolutions gives evidence for an improved representation of local precipitation processes at a horizontal resolution of 50 m compared to 450 m. Additionally, differences of precipitation between 2830 m above sea level and the ground indicate that near-surface processes are active in the model.

1 Introduction

In many regions of the world, e.g., the Alps or the Californian Sierra Nevada, snow is the main source of fresh water. Additionally, it is an important resource for hydropower and is crucial for winter tourism in skiing areas . Thus, especially in a changing climate, it is essential to improve the understanding of processes forming seasonal snow cover. Improving the ability of weather forecast models to represent the spatial variability in snowfall is further crucial to efficiently manage fresh water and hydropower. Moreover, as snow is a potential danger in terms of avalanches, improved knowledge about the distribution of snow is crucial for avalanche forecasting and prevention.

Snow accumulation patterns at a mountain-range scale are known to be strongly dependent on blocking and lifting processes including large-scale orographic precipitation enhancement (e.g., Houze Jr.2012; Stoelinga et al.2013), which is related to large-scale atmospheric circulation. However, for a long time little knowledge was available about the spatial distribution of snow on a mountain-slope or river-catchment scale. Only in recent years have improvements in technology allowed the investigation of mountain-slope scale snow distribution (e.g., Deems et al.2006; Prokop2008; Grünewald et al.2010). Terrestrial and airborne laser scanning reveal annually persistent patterns of peak-of-winter snow accumulation distribution on river-catchment scales , which is found to be consistent with few dominant snowfall events of the season. Reported scale breaks in fractal analysis of snow accumulation patterns are mainly at scales of <100 m and represent the occurrence of a change in dominant processes (e.g., Deems et al.2008). On very small scales snow accumulation patterns are assigned to wind redistribution of snow (e.g., Mott et al.2011; Vionnet et al.2017). Vegetation effects were found to be dominant at small scales and terrain effects dominate on scales of up to 1 km . Different dominant scales are reported for different slope expositions relative to the wind direction . Furthermore, could show that snow accumulation smooths the underlying terrain, reducing the small-scale spatial variability in topography. While most studies addressed variability in snow accumulation, the combined scale analysis of snow accumulation and snow precipitation patterns by reveals much smoother patterns in snow precipitation at about 300–600 m above ground compared to final snow accumulation at the ground on scales of up to 2 km. This stresses the importance of pre-depositional near-surface and post-depositional processes for snow accumulation patterns.

Driving processes of snow accumulation on the mountain-ridge scale were addressed in numerous studies, which reveal two main pre-depositional processes. On the one hand, mountain-ridge-scale precipitation and accumulation are influenced by local cloud dynamical processes . On the other hand, particle–flow interactions (i.e., the influence of the local flow field on the pathways of snow particles and particle distribution in the air) determine snow accumulation patterns in mountainous terrain . On the mountain-ridge scale, documented the occurrence of a local event of orographic snowfall enhancement. In their case study, the presence of a low-level cloud gives evidence for precipitation enhancement favored by the seeder–feeder mechanism (e.g., Bergeron1965; Purdy et al.2005). On similar scales, preferential deposition was found to cause enhanced snow accumulation on leeward slopes (e.g., Mott et al.2010; Mott and Lehning2010). However, snow depth measurements in very steep terrain and corresponding local flow field measurements reveal even more complex particle–flow interactions than previously suggested by model studies. On even smaller scales the main driver of snow accumulation patterns is post-depositional snow transport by drifting and blowing snow, which is dependent on local topographic features and wind gusts .

Complex terrain–flow–precipitation interactions (i.e., the effect of terrain-induced flow field variations on precipitation formation and distribution), especially on the mountain-ridge scale, still leave the relative importance of the different pre-depositional processes for snow accumulation and the frequency of occurrence barely known . Running a coupled simulation of the snowpack model Crocus and the atmospheric model Meso-NH in large-eddy simulation (LES) mode, addressed the question of the relative importance of these different processes including snow redistribution by wind. Their results show that post-depositional snow transport dominates snow accumulation variability but leaves the question of the relative importance of pre-depositional processes open.

Given the small scale of these processes, their relative importance may be addressed either based on very-high-resolution numerical simulations or based on spatially highly resolved precipitation measurements. Therefore, accurate model results and radar measurements at high resolution are essential. Both, however, are challenging to achieve and very-high-resolution simulations are still rare, especially over complex terrain. Remote-sensing techniques are the most important methods to obtain high-resolution spatial measurements of atmospheric properties at different atmospheric levels. They permit us to gain information about both the small- and the large-scale properties of the atmospheric processes. Of particular importance is weather radar due to its wide coverage, fine spatial resolution and interaction of microwaves with the precipitation. These properties have been used to infer orographic precipitation enhancement, particularly in the case of liquid precipitation .

In this study, we present very-high-resolution WRF simulations, which are forced by 2.2 km resolution Consortium for Small-scale Modeling (COSMO) analysis and high-resolution radar estimates making use of the recently renewed MeteoSwiss radar network and its adequate technical performances, which allow the observation of precipitation in a challenging, complex alpine environment. Combining the COSMO–WRF simulations with operational radar measurements, we perform a variability analysis for snow precipitation at a regional to mountain-ridge scale to address the following question: to what degree is snow precipitation variability represented by very-high-resolution WRF simulations?

2 Data and methods

## 2.1 WRF model setup

Atmospheric simulations are performed with the non-hydrostatic and fully compressible Weather Research and Forecasting (WRF) model version 3.7.1 for the region of eastern Switzerland (Fig. 1). Simulations are set up with four one-way nested domains (d01–d04, Fig. 1). Domain d01 has a horizontal resolution of 1350 m with 40 vertical levels and covers a region of about 250 km × 320 km including eastern Switzerland and a portion of the neighboring countries (Fig. 1, Table 1, Supplement S1). The three nests have horizontal resolutions of 450, 150 and 50 m using a nesting ratio (dxparent∕dxnest) of 3. Domains d02–d04 have 40, 60 and 90 vertical levels, respectively, with the model top at 150 hPa using a preliminary version of vertical nesting . A total of 20 and 40 vertical levels refine the whole atmosphere in domains d03 and d04, respectively. A total of 10 vertical levels in d04 are introduced to additionally refine the boundary layer. To make sure that there is plenty of domain for the model to adapt to the refined topography, domain d02 is shifted toward the eastern boundary of domain d01 as dominant wind directions are from the northwesterly and southerly directions. Domain d02 covers the central northern part of the Grisons, while domains d03 and d04 cover the surroundings of Davos and the upper Dischma valley, respectively (Fig. 1). Simulations are performed for three snow precipitation events on 31 January, 4 February and 5 March 2016 (Sect. 2.5).

Figure 1Overview over the study area in the eastern part of Switzerland surrounding Davos. WRF simulation domains (d01–d04, dark red to yellow) and evaluation domains (blue) give information on the simulation and evaluation setup. The 18 meteorological stations (red triangles) are within or very close to the regional domain. The two stations Dischma Moraine and FLU2, which are used to validate the model, are within the domain Dischma. The operational weather radar is located on the Weissfluh summit at approximately 2830 m above sea level (m a.s.l., blue pentagon). Coordinates in the right panel are in Swiss coordinates CH1903LV03 (unit: m). Shaded topography: dhm25 ©2018 swisstopo (5 740 000 000).

The parent domain is run with a planetary boundary layer (PBL) scheme (non-large eddy simulation (non-LES) mode), while the three nests are run in the LES mode. No strong differences were found when running domains d02 and d03 in non-LES mode (not shown). Therefore, and as we are interested in having an as good as possible representation of small-scale winds, we decided to run our simulations in the LES mode for all nested domains. Domains d02 and d03 are within the “gray zone” (Wyngaard2004). There are approaches omitting simulations in the gray zone by the choice of a higher grid refinement ratio , which would be worth a sensitivity study. However, we use the well-tested 1 : 3 grid refinement ratio and keep our model setup consistent with the very-high-resolution simulations by , except that they perform separate simulations for the non-LES and LES domains, while we run a nested simulation with one-way feedback for all four domains. Running a nested simulation of the non-LES and LES domains turned out to be necessary for precipitation to evolve properly in the LES domains, as hydrometeors cannot be used as a boundary condition for the parent domain but are fed to nested domains in WRF simulations. Sub-grid-scale turbulence is parametrized by the 1.5 order turbulent kinetic energy closure . For the non-LES setup the Yonsei University PBL parameterization (YSU PBL; Hong et al.2006), which is considered to be one of the schemes showing the best performance over complex terrain , is used. An adapted version of YSU PBL was shown to perform even better when taking into account sub-grid-scale variability in the terrain . However, given our high model resolution we decided to keep the model simple and run the simulations with the standard YSU PBL. Land use data are taken from the CORINE dataset and translated to the U.S. Geological Survey conventions . Soil type is set to silty clay loam for the whole domain. The link between the soil, which is modeled by the Noah land-surface model with multi-parameterization options (Noah-MP, Niu et al.2011; Yang et al.2011), and the atmosphere is given by the MM5 Monin–Obukhov surface layer model , which is based on the Monin–Obukhov similarity theory (Obukhov1971). For microphysics the Morrison two-moment precipitation scheme , which was found to be one of the schemes that most adequately simulates snow precipitation over complex terrain , is used. Details about processes in the Morrison parameterization are given in Appendix Appendix A. An investigation of different microphysical parameterizations would be interesting but is beyond the scope of this study. Given the high horizontal resolution, no sub-grid parameterizations for cumulus clouds are used.

The 2.2 km horizontally resolved Consortium for Small-scale Modeling (COSMO-2) analysis by MeteoSwiss is used as initial and boundary conditions for the parent domain. For COSMO-2 analysis data to be readable by the WRF preprocessing system a regridding of the rotated COSMO coordinates to latitude–longitude coordinates is required. COSMO preprocessing, model adaptations and details about the model simulations are given in .

Table 1Setup for the four nested domains (d01–d04) used in the WRF simulations. For the planetary boundary layer (PBL) the simulation mode is given, distinguishing between non-large eddy simulation (LES) and LES settings. For non-LES settings the PBL scheme is given. Additionally, the sub-grid-scale (SGS) turbulence parameterization is given for all four domains. dx and dy give the horizontal resolution. Vertical levels gives the number of vertical levels in the different domains. Vertical levels (<1000 m) gives the number of vertical levels in the lowest 1000 m of the atmosphere. The lowest 21 model levels for all four domains are given in Table S1. The time step (dt) and the maximum slope angle are given for simulations with four (14) smoothing cycles.

1 YSU: Yonsei University PBL scheme. 2 TKE: turbulent kinetic energy.

Topography in the model is based on the ASTER Global Digital Elevation Model V002 with a resolution of 1 arcsec . Terrain smoothing has been applied for all domains due to the very steep terrain in the simulation area. Four cycles of the WRF 1–2–1 smoothing (i.e., a moving window filter with a window length of 3 and weights of 1 : 2 : 1 for the grid points i−1, i and i+1) are applied to all four domains to keep all slopes in domain d04 (50 m horizontal grid spacing) below 45. Additionally, the boundaries of the parent domain are smoothed to match COSMO topography . Test simulations are run with 14 cycles of WRF 1–2–1 smoothing, which allows for a longer computational time step and therefore saves computational time (Table 1). Maximum slope angles for all domains and different smoothing are given in Table 1. Simulations with different precision of topography further allow us to address the importance of the representation of topography in the model. To allow the simulations to adapt to higher-resolution topography, domains d01–d04 are run with a spin-up of 43, 19, 7 and 1 h, respectively.

As the snow cover in complex alpine terrain is likely rougher than for a flat field and to account for non-resolved topography and additional smoothing, snow surface roughness length has been changed to 0.2 m. The chosen roughness length is much larger than roughness lengths assumed by , for example. However, grid spacing in our simulations is larger and the roughness length is chosen such that it accounts for roughness elements in complex terrain (e.g., large rocks) and non-resolved topography, which are assumed to have an average size of about 2 m. This estimate is based on a comparison of a 2 m digital terrain model (DTM-AV © 2018 swisstopo (5 704 000 000)) to a 25 m resolution digital elevation model (dhm25 © 2018 swisstopo (5 740 000 000)), which reveals an average difference on the order of 2.5 m for bare ground conditions in domain d04 between 2200 and 2700 m a.s.l. Hence, the estimate of 2 m is rather conservative but takes into account smoothing of the terrain by the snow cover.

For the model validation (Sect. 3.1), WRF variables, using model output of the four grid points surrounding the station, are linearly interpolated to the coordinates of the meteorological station (see Sect. 2.2). Alternatively, the eight neighbors of the grid point closest to the station could be used . Temperature is corrected for elevation due to terrain smoothing using a moist-adiabatic temperature gradient of −0.0065 km−1. Modeled wind speeds are extrapolated to the measurement height by applying a logarithmic wind profile, as wind measurements at the automatic weather stations are not taken at 10 m but 4 or 5 m above ground (Sect. 2.2). This is a rough approximation given the assumption of a neutral atmosphere. For simulation domains d01–d03 10 m wind speeds are extrapolated to the elevation of the sensor above the snow cover, while for domain d04 wind speeds at the lowest model level (approximately 3 m above ground) are used for the extrapolation. The dynamic reference roughness length is chosen to be 0.2 m (corresponding to the surface roughness length in the model simulations). For wind direction comparisons wind directions at 10 and 3 m above ground are chosen for the simulation domains d01–d03 and d04, respectively. As a reference COSMO-2 variables of the closest grid point to the station are included in the model validation and hence in Figs. 2 and 3. The 2 m temperature and 10 m wind speed of COSMO-2 are corrected for elevation using the same procedure as for the WRF simulations.

## 2.2 Automatic weather stations

Snow depth measurements from a total of 18 automatic weather stations in the central northern part of the Grisons (Fig. 1) are used. Two stations (Dischma Moraine and Dischma Dürrboden) were installed as part of the Dischma Experiment (DISCHMEX), in which processes of snow accumulation and ablation in the Dischma valley near Davos (Switzerland) are addressed . A total of 16 stations are part of the Intercantonal Measurement and Information System (IMIS). The 18 stations are located between 1560 and 2725 m a.s.l. The stations measure snow depth in addition to the standard meteorological parameters. All stations have shielded temperature and humidity sensors but are unheated. Biased temperatures around midday and occasional data gaps due to iced instruments may therefore occur . Two stations (Dischma Moraine and FLU2), which are located in the WRF domain d04 with a horizontal grid spacing of 50 m, are used for the model validation. The variables evaluated are 2 m temperature, 2 m relative humidity, wind speed and wind direction. Wind measurements at IMIS stations are taken about 5 m above ground, while the wind sensor at station Dischma Moraine is located at about 4 m above ground.

## 2.3 Operational weather radar data

Weather radar datasets employed in the presented analyses are acquired by the MeteoSwiss operational radar located at the Weissfluh summit (2850 m a.s.l.), in the proximity of Davos. It is a dual-polarization Doppler weather radar, providing complementary information about the detected hydrometeors by considering their interaction with the incident electromagnetic radiation in both horizontal and vertical polarization planes. This complementary information leads to an enhanced clutter detection, which makes the radar measurements in such a complex mountainous terrain significantly more reliable. The polarimetry also makes it possible to identify the type of hydrometeors , which allows us to be confident that in the zone of interest for the presented study we deal with solid precipitation, consisting mostly of aggregates and crystals and partly of rimed ice particles.

The radar operates in 5 min cycles during which it scans the surrounding atmosphere by performing complete rotations at 20 different elevations, from −0.2 to 40 . Operationally, the size of a radar sampling volume is 500 m in range, whereas the size observed in the perpendicular plane depends on the half-power beamwidth and increases with range. The acquired data undergo an elaborated procedure of corrections . Before the quantity of precipitation at the ground level is estimated by averaging over 1 km2, the observations are corrected for the vertical profile of reflectivity (VPR) with the weight assigned to volumes being inversely proportional to their height above the ground .

In the framework of our study, rather than relying on the operational radar product, we use data with the highest available resolution of 83 m in range. We also adopted a more conservative, nonoperational method of clutter identification, which relies exclusively on the polarimetry and leaves very little residual clutter, however, sometimes at the expense of removing some precipitation. Given that we consider only radar measurements at low elevation angles in the vicinity of the radar and that the bright band is not present in our case studies (all radar measurements are from above 2800 m a.s.l. during the winter season), the observations are not corrected for the VPR. Furthermore, given the strong influence of wind on the snow precipitation, we restrict our precipitation estimate to only four elevations, from the second to fifth (0.4, 1, 1.6, 2.5), avoiding the first one, judged to contain too little information due to the abundant rejected ground clutter areas.

Polarimetry helps to identify non-meteorological scatterers, to distinguish between different types of hydrometeors, to correct for signal attenuation and to make quantitative estimates of intense to heavy rainfall. For snowfall measurements it is common to use reflectivity Z at horizontal polarization and convert it into snow water equivalent S using a so-called ZS relationship :

$\begin{array}{}\text{(1)}& Z=\mathrm{100}{S}^{\mathrm{2}}.\end{array}$

The coefficients used in this formula account for the dielectric properties and fall velocities of snow and convert reflectivity Z in snow water equivalent S. The radar provides an indirect estimate of snowfall rather than a direct measurement. Applied to each radar sampling volume scanned by the four selected elevations in the zone of interest (up to 40 km around the radar), the formula gives an estimate of liquid precipitation equivalent in the three-dimensional volume. By vertically averaging estimates from the four elevation sweeps using equal weights, we obtain the estimate of precipitation in polar (range, azimuth) coordinates at a flat plane at the height level of the radar. These estimates are summed up over 24 h to obtain the accumulation maps used in the study.

Further on, the polar accumulation maps are resampled by means of the bilinear interpolation to the Cartesian grid of the regional domain (450 m resolution) and the local domain (300 m resolution). The obtained Cartesian maps are finally processed to remove the residual clutter using a 3×3 median filter, partly or entirely. The former means that only the isolated high values in the original map are replaced with the corresponding value of the filtered map, at the positions where the difference between the original and the filtered map appears to be larger than 5 mm (hereafter “partly filtered”). The latter means that the entire map is influenced by the median filtering (hereafter denoted as “entirely filtered”).

## 2.4 Autocorrelation and variogram analysis

To investigate the variability in snow precipitation and accumulation patterns and their relation to topography, a scale analysis, based on two-dimensional (2-D) autocorrelation maps and variograms, is performed. Two-dimensional autocorrelation maps and variograms are further used to relate variability in radar and WRF precipitation. Given the resolution restriction by the radar measurements (Sect. 2.3), we analyze two different domains using horizontal resolutions of 450 and 300 m, respectively. The domain with a resolution of 450 m covers an area of about 58 km × 56 km centered over the radar on the Weissfluh summit (hereafter regional domain, Fig. 1). The domain with a resolution of 300 m covers an area of 24 km × 21 km to the south of Davos (Switzerland) including the Dischma valley (hereafter local domain). Radar data (300 m resolution) and WRF precipitation at three resolutions (450, 150 and 50 m) are additionally evaluated for domain Dischma to address the influence on the spatial resolution of variability. Domain Dischma covers the upper Dischma valley with an extent of 9 km × 9 km.

To produce variograms the semivariance ($\stackrel{\mathrm{^}}{\mathit{\gamma }}$) is calculated at 50 logarithmic lag distance bins (h, i.e., a set of distance ranges) by

$\begin{array}{}\text{(2)}& \mathit{\gamma }\stackrel{\mathrm{^}}{\left(}h\right)=\frac{\mathrm{1}}{\mathrm{2}|N\left(h\right)|}\sum _{\left(i,j\right)\in S\left(h\right)}\left({a}_{j}-{a}_{i}{\right)}^{\mathrm{2}},\end{array}$

where S(h) are the point pairs (i,j) and N(h) gives the number of point pairs of the evaluated variable a. WRF and radar snow precipitation and topography are evaluated at 450 and 300 m resolutions with maximum lag distances of 25 and 10 km, respectively. Variograms for domain Dischma are calculated with a maximum lag distance of 5 km. Minimum numbers of point pairs in one lag distance bin for the local and regional domain are 18 317 and 8035, respectively. For domain Dischma the minimum number of point pairs is between 677 and 55 419, depending on the resolution.

To determine scaling properties, an empirical log-linear model is fit to the variogram by least-square optimization . The model used is not a valid variogram model but used to describe the experimental variograms and chosen to be consistent with Schirmer et al. (e.g., 2011). For all variograms three empirical log-linear models are fit:

$\begin{array}{}\text{(3)}& y\left(x\right)=\left\{\begin{array}{ll}{\mathit{\alpha }}_{\mathrm{1}}×\mathrm{log}\left(h\right)+{\mathit{\beta }}_{\mathrm{1}},& \text{for}\phantom{\rule{0.25em}{0ex}}\mathrm{log}\left(h\right)<{l}_{\mathrm{1}};\\ {\mathit{\alpha }}_{\mathrm{2}}×\mathrm{log}\left(h\right)+{\mathit{\beta }}_{\mathrm{2}},& \text{for}\phantom{\rule{0.25em}{0ex}}{l}_{\mathrm{1}}\ge \mathrm{log}\left(h\right)<{l}_{\mathrm{2}};\\ {\mathit{\alpha }}_{\mathrm{3}}×\mathrm{log}\left(h\right)+{\mathit{\beta }}_{\mathrm{3}},& \text{for}\phantom{\rule{0.25em}{0ex}}\mathrm{log}\left(h\right)\ge {l}_{\mathrm{2}};\end{array}\right\\end{array}$

using the constraint that each log-linear model needs to contain a minimum of four data points and the continuity constraint(s)

$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{1}}\mathrm{log}\left(l\mathrm{1}\right)+{\mathit{\beta }}_{\mathrm{1}}& ={\mathit{\alpha }}_{\mathrm{2}}\mathrm{log}\left(l\mathrm{1}\right)+{\mathit{\beta }}_{\mathrm{2}},\\ {\mathit{\alpha }}_{\mathrm{2}}\mathrm{log}\left(l\mathrm{2}\right)+{\mathit{\beta }}_{\mathrm{2}}& ={\mathit{\alpha }}_{\mathrm{3}}\mathrm{log}\left(l\mathrm{2}\right)+{\mathit{\beta }}_{\mathrm{3}}\end{array},\end{array}$

where ${\mathit{\alpha }}_{\mathrm{1},\mathrm{2},\mathrm{3}}$ and ${\mathit{\beta }}_{\mathrm{1},\mathrm{2},\mathrm{3}}$ are the slopes and intercepts of the three log-linear models, respectively. Scale breaks (l1, l2) are the lag distances of the intersections of the first and second and second and third log-linear models, respectively. Scale breaks were previously found to determine the scale of a change of dominant processes (e.g., Deems et al.2006). To address the variability with respect to the overall variability in the respective domain, all variograms are normalized by the total domain-wide variance.

Two-dimensional autocorrelation is calculated based on Pearson's correlation coefficient r of all grid point pairs for a maximum lag distance of ±40 grid points in the x and y directions. This results in maximum lag distances of 18 km for the regional domain.

## 2.5 Snowfall events

This study is based on three precipitation events in winter 2016. On 31 January 2016 the Azores high and a low-pressure area over Scandinavia induced westerly flow over central Europe and relatively mild temperatures with about −3C at 2500 m a.s.l. A shift of the Azores high toward northern Spain and a trough over eastern Europe led to a change in wind direction toward northerly advection and a decrease in temperature (about −12C at 2500 m a.s.l.) on 4 February 2016. On 5 March 2016 a low-pressure area over France, which is part of a large depression area over central Europe, causes southerly advection over Switzerland. Temperatures were about −7C at 2500 m a.s.l. Given the relatively high temperatures on 31 January 2016, which resulted in quite substantial liquid precipitation at the lowest elevations, total (solid and liquid) ground precipitation is evaluated. This does not make a big difference for the precipitation events on 4 February and 5 March 2016 but is essential for the precipitation event on 31 January 2016. For precipitation patterns at the elevation of the radar (2830 m a.s.l.) we only analyze solid precipitation from WRF.

Figure 2Comparison of (a–c) 2 m temperature (C), (d–f) 2 m relative humidity (%), 4 m (g–j) wind direction () and (k–m) wind speed (m s−1) at the station Dischma Moraine (black) to WRF simulations interpolated to the coordinates of the station Dischma Moraine for the three precipitation events on 31 January, 4 February and 5 March 2016 for all four simulation domains (d01: dark green, d02: light green, d03: light blue, d04: dark blue). For comparison COSMO-2 is added for the closest grid point (dashed gray). The 2 m temperature in WRF and COSMO is corrected for elevation based on a moist-adiabatic temperature gradient.

Figure 3As in Fig. 2 but for the station FLU2 on the Flüela Pass with 5 m wind speed and direction.

3 Results and discussion

## 3.1 Point validation of WRF simulations

2 m air temperature and relative humidity, and 4 or 5 m wind speed and direction at two stations (Sect. 2.2) are compared to WRF to validate the model (Figs. 2 and 3).

For both stations 2 m temperature matches reasonably with observations, although especially for the precipitation event on 4 February 2016 substantial temperature deviations occur around midday (Figs. 2a–c and 3a–c). Deviations of the WRF model from station measurements during midday are likely caused by errors in station measurements due to radiative heating of the multiplate shielded temperature sensors (Huwald et al.2009, Sect. 2.2).

Relative humidity shows partially good agreement but shows a strong temporal variability (Figs. 2d–f and 3d–f). WRF is generally able to capture the main drops in relative humidity at the two investigated stations but it introduces additional drops compared to measurements. The microphysics parameterization is originally developed for simulations with a coarser resolution, which produce fewer vertical motions. Thus, the introduction of a higher variability in relative humidity in our WRF simulations may be due to strong subsidence and lifting, which lead to an overestimation of adiabatic cooling or warming and hence to an overestimation of humidity generation or decay. Additionally, differences between modeled and measured relative humidity may be due to measurement uncertainties.

Simulated wind direction shows good agreement with measured wind direction (Figs. 2g–i and 3g–i). In complex terrain, simulations are often limited to resolutions, which are too coarse to resolve smaller-scale terrain features that affect near-surface wind direction (e.g., due to a lack of high-resolution terrain data, or computational resources), and thus cannot accurately capture changes in wind direction close to the surface where weather stations are located. Good agreement in wind direction modeling in our COSMO–WRF simulations in complex terrain is likely due to the high resolution of topography. For some cases wind directions in the WRF simulations additionally improve for higher resolutions, although for others terrain smoothing is likely to have adverse effects on modeled wind direction.

Compared to the good agreement of wind direction, wind speeds show only partially good agreement with station measurements (Figs. 2k–m and 3k–m). Wind speeds were found to strongly depend on the sub-grid-scale turbulence parameterization and a strong overestimation of wind speeds was observed for different simulation setups (not shown). Applying the improved nonlinear sub-grid-scale turbulence parameterizations leads to instabilities in the current model setup. The use of a snow surface roughness length of 0.2 m, representing the combined roughness of snow and surface features (e.g., boulders and rocky outcroppings; Sect. 2.1), compared to simulations with a standard WRF roughness length of snow of 0.002 m, could partially reduce overestimated wind speeds . While we address non-resolved topography based on an increased snow surface roughness length, another approach to improve wind speeds in WRF simulations has been introduced by , who use a sink term in the momentum equation based on sub-grid-scale topography. They demonstrate the ability of their approach to improve surface wind speeds. However, with this approach, the effect of the sub-grid-scale topography is only included for simulations using a PBL parameterization. As in our model setup a PBL parameterization is only applied for domain d01; we address the non-resolved topography by increasing the surface roughness, which allows us to include the effect of non-resolved topography for all four simulation domains. In addition, our simulations are run over snow-covered terrain, which implies that the standard roughness length for snow used in WRF (on the order of 10−3 m) is 2 orders of magnitude lower than roughness lengths representative of the scale of complex terrain in our simulations (on the order of 10−1 m). Applying the PBL version of might be a possibility to reduce excess wind speeds in domain d01, which might also impact wind speeds in the domains d02–d04. However, such a sensitivity study is out of scope of the present investigation.

Based on our approach COSMO–WRF still simulates excess wind speeds for the two precipitation events on 31 January and 4 February 2016 and for the precipitation event on 5 March 2016 at station FLU2 for a few hours around noon. The overestimation on 31 January and 4 February 2016 is assumed to be connected to the upwind location of both stations during these two precipitation events, as speedup over windward slopes and ridges is a known problem . Hence, the exact location of the station relative to the ridge is important to verify wind speeds. Furthermore, local terrain features upstream of the station may disturb the wind field. For example, the station Dischma Moraine is located on a moraine on the northern side of the ridge between Piz Grialetsch and Scalettahorn. The station FLU2 is located on the northern side of Flüela Pass above a small rock face and to the east of a terrain knoll. Such terrain features, while not represented in the model, may strongly reduce wind speeds in reality. Another potential cause of disagreement is that station measurements are prone to measurement uncertainties, and riming of the unheated instruments may lead to an underestimation of wind speeds .

Generally, reasons for an overestimation of wind speeds may be manifold. An exact estimation of wind speeds at stations in the model is not expected due to unresolved topographical features in the complex terrain of our study site. An additional source of uncertainty – though unlikely to be on the order of the strong excess wind speeds – is the extrapolation of wind speeds based on the assumption of a neutral atmosphere. While different potential causes of wind speed overestimation are discussed above, actual reasons for deviations in wind speed remain unknown.

While the model is designed such that it develops independently (given its fetch distances and spin-up times, Sect. 2.1), a poor representation of the large-scale gradients in the COSMO-2 input might not be corrected by COSMO–WRF. The investigated variables do not show a consistent signal of improvement nor a consistent signal of worsening with higher resolutions and with respect to the COSMO-2 input. Similarly but not consistently in phase with COSMO–WRF, COSMO-2 shows a good agreement with station measurements for certain cases, while it shows worse performance for other cases. Given these inconsistencies between cases and for COSMO-2 input, a poor representation of the large-scale gradient in COSMO-2 is an unlikely reason for bad model performance.

Overall, we show that the presented simulation setup captures temperature, relative humidity and wind conditions in complex terrain at two stations by a certain degree. Temperature deviations around midday are likely due to measurement uncertainties. Wind speeds tend to be overestimated, especially on the windward side of the mountain ridges. This shows that even for very high model resolutions the point performance of the model with respect to wind speeds remains challenging. However, the model shows a good performance in the representation of the local wind directions, which is likely a consequence of the improved representation of topography.

## 3.2 Spatial snow precipitation and accumulation patterns

Radar precipitation maps of the regional domain covering an area of about 58 km × 56 km centered over the radar on the Weissfluh summit (Fig. 1) tend to show precipitation patterns dependent on wind direction (Figs. 2g–i, 3g–i and 4d–f) (Fig. 4). The precipitation field on 31 January 2016 shows a strong south–north gradient (Fig. 4a), while the precipitation field on 4 February 2016 shows a more homogeneous distribution (Fig. 4b). For the precipitation event on 5 March 2016 radar precipitation maxima are observed over the mountain ridges in the southern part of the domain (Fig. 4c). Although our regional domain does not represent a cross section across the whole alpine mountain range, a north–south (south–north) precipitation gradient for southerly (northerly) advection is apparent. This is in good agreement with large-scale orographic precipitation enhancement , which favors precipitation on the upwind side of a mountain range due to topographically induced lifting and a drying due to sinking air masses downwind of the mountain range.

These large-scale patterns of orographic precipitation enhancement are partially captured in the WRF simulations (Fig. 4d–f). In particular, for southerly advection (precipitation event on 5 March 2016) this large-scale effect is well represented in COSMO–WRF, in which precipitation maxima occur over mountain ridges in the southern part of the domain and a north–south precipitation gradient is present. For northerly to northwesterly advection (precipitation events on 31 January and 4 February 2016), however, snow precipitation maxima in the WRF simulations are shifted eastward compared to radar precipitation estimates, i.e., toward the outflow boundary.

Microphysics and precipitation dynamics in the model are likely to be a limiting factor in terms of small-scale precipitation patterns. Disagreement between radar and WRF precipitation patterns may further be connected to the strong terrain smoothing in the model. Despite the high resolution of our simulations, slope angles are relatively low, with maximum slope angles of 35.2 in the regional domain due to the application of terrain smoothing (Table 1). Given even lower slope angles in domain d01, precipitation fed to domain d02 may already be too weak and thus needs to develop within domain d02. As mountains in the northwestern part of the domain are shallower than mountains in the southeastern area (Fig. 1), lifting condensation may not be strong enough in the northwestern area of the domain, leading to precipitation generation further downstream in the domain, where steeper and higher mountains may even lead to too strong precipitation enhancement. Additionally, if the tendency of overestimated wind speeds sustains up to higher atmospheric levels in the model, this may lead to an overestimation of the advection of hydrometeors in the microphysics scheme . This would result in a downstream shift of the precipitation maximum. However, we do not expect this to have a strong impact on the regional-scale precipitation distribution. Thus, there are likely additional reasons for the observed downstream shift of WRF precipitation compared to radar precipitation, which remains difficult to explain.

Figure 5Domain-wide 24 h precipitation statistics for the regional domain (450 m resolution, Fig. 1) for the three precipitation events on 31 January, 4 February and 5 March 2016. Gray colors show entirely filtered radar precipitation. WRF precipitation at 2830 m above sea level (m a.s.l.) for simulations with weak terrain smoothing (Sect. 2.1) and strong terrain smoothing are given in blue and violet, respectively. Orange (red) shows box plots of WRF total ground precipitation for weak (strong) terrain smoothing. Radar precipitation and WRF precipitation at 2830 m a.s.l. are masked (as shown in Fig. 4).

New snow depth measured at 18 automatic weather stations in the regional domain (Fig. 4k–m) over 24 h shows a distinct elevation gradient, which is quite well represented by WRF total ground precipitation (Fig. 4g–i). For 31 January and 5 March 2016 the large-scale precipitation trend observed in the radar data is generally represented in station measurements. On 4 February 2016 station measurements suggest a precipitation peak in the upper Dischma valley (lower left quadrant in Fig. 4l), which agrees with WRF simulations. Radar estimates, however, show a more homogeneous distribution of precipitation on 4 February 2016. Snow depth changes at the stations are very local and strongly affected by wind redistribution of snow, which may disturb the large-scale gradient. Additionally, the distribution of stations is not homogeneous over the regional domain and fewer stations are available in the western part of the domain.

The visual comparison of radar and WRF precipitation patterns for all three events (Fig. 4) reveals that precipitation patterns are influenced by wind direction and topography. Large-scale precipitation patterns are in agreement with station measurements, although the latter are strongly influenced by the local wind field and snow redistribution processes.

## 3.3 Mean variability

Radar precipitation distributions at 2830 m a.s.l. in the regional domain (450 m resolution, Fig. 5) show a larger interquartile range (IQR) than radar precipitation in the local domain (300 m resolution, Supplement S2), confirming that local precipitation is more uniform than regional precipitation. Radar median precipitation over 24 h is on the order of 10 to 20 mm water equivalent for all three precipitation events in the regional domain. The median of radar precipitation in the local domain can be both higher or lower than in the regional domain.

Although radar estimates are based on a reference SZ relationship, the employed formula (Eq. 1) is not immune to potential estimation errors. Therefore, despite reasonably assuming that the potential estimation errors should not significantly influence the variability in and the relative intensity of the precipitation fields, we consider potential inaccuracies in our interpretations.

For the precipitation events on 31 January and 4 February 2016 median precipitation at 2830 m a.s.l. in the COSMO–WRF simulations is in reasonable agreement with radar median precipitation (Fig. 5), even though WRF and radar precipitation patterns are different (Sect. 3.2). However, for the precipitation event on 5 March 2016 the median of precipitation in the regional domain is higher in WRF simulations compared to radar measurements, while the large-scale precipitation gradient is in good agreement (Fig. 4c and f). The IQR of WRF precipitation is generally larger compared to the IQR of radar precipitation and the domain-wide WRF precipitation distribution has longer tails compared to the radar precipitation distribution. In the local domain these tendencies are preserved (Supplement S2). Median precipitation is higher in the WRF simulations compared to radar estimates on 31 January 2016 in addition to 5 March 2016, for which the deviations increase. This supports the hypothesis that the model tends to overestimate precipitation for higher resolutions with steeper and more complex topography.

The domain-wide median and IQR of precipitation at 2830 m a.s.l. in WRF simulations with weaker and stronger terrain smoothing (Sect. 2.1) are similar with a slight tendency of higher median values for weaker smoothing, indicating that the accuracy of topography does not have a strong influence on the domain-wide statistics of precipitation on the regional scale. Enhanced precipitation for weaker terrain smoothing compared to stronger terrain smoothing could be explained by enhanced precipitation production due to steeper topography.

An overestimation of precipitation in WRF simulations was previously reported and could be due to various reasons. and among others report that WRF tends to show a stronger overestimation of precipitation for higher model resolutions compared to coarser model resolutions. Additionally, they document a dependency on the intensity of precipitation. An overestimation of orographic precipitation enhancement in more complex terrain or an overestimation of moisture in the model were further reported by . An overestimation of orographic precipitation enhancement would be in agreement with a stronger overestimation of precipitation for the local domain compared to the regional domain and for weaker smoothing compared to stronger smoothing. Furthermore, it is likely to occur for simulations with high horizontal resolution as higher peaks and steeper slopes are preserved . Compared to a shallow topography, higher peaks and steeper slopes may cause stronger lifting and subsidence, which is also a likely cause for additional drops in relative humidity in WRF compared to measured relative humidity (Sect. 3.1, Figs. 2 and 3). This tendency seems to only apply for the highest elevations. For lower elevations strong smoothing may result in elevation differences, which are too small for precipitation to evolve by lifting condensation (Sect. 3.2). As additional reasons for precipitation overestimation in WRF, an overestimation of precipitation in the driving model and underlying land use characteristics were mentioned. The latter was, however, previously found to only have a weak influence on the precipitation amount (Pohl2011). Humidity in COSMO-2 is an unlikely reason as COSMO-2 shows a tendency of underestimating relative humidity compared to station measurements (Figs. 2 and 3). Even though there are many possible reasons for an overestimation of precipitation in WRF, the estimation of solid precipitation from radar measurements is also subject to uncertainties (e.g., Cooper et al.2017). Given uncertainties in radar precipitation estimates, the comparison of median domain-wide precipitation should be taken with care. An in-depth analysis of spatial variabilities is given in Sect. 3.4.

At the ground level WRF precipitation tends to show higher median values of precipitation compared to WRF precipitation at 2830 m a.s.l. for both domains. The IQR is similar. From this we hypothesize that there are precipitation formation or enhancement processes taking place between the elevation of the radar and the ground. This is in good agreement with the fact that near-surface processes can strongly enhance snow precipitation (e.g., riming). Overall, this analysis shows that WRF tends to overestimate domain-wide precipitation and precipitation variability at 2830 m a.s.l. compared to radar estimates.

Table 2Large-scale linear trends of entirely filtered radar and WRF precipitation patterns on the regional domain (450 m horizontal grid spacing, Fig. 1). WRF precip. at 2830 m a.s.l. refers to solid precipitation in WRF simulations at 2830 m above sea level and WRF total ground precip. refers to the total (solid and liquid) precipitation at the ground level. Orientation gives the direction of the slope and intensity the strength of inclination. A value of 0 would indicate a slope pointing toward the east. WRF snow precipitation is from simulations with weak terrain smoothing (Sect. 2.1).

## 3.4 Spatial variability

To address spatial patterns and variability in precipitation, a scale analysis is performed augmented with a 2-D autocorrelation analysis (Sect. 2.4). Given the overestimation of precipitation in the model and the large differences in domain-wide variability between the model and radar precipitation estimates (Sect. 3.3), all variograms are normalized by the domain-wide variability, which allows analysis of spatial patterns with respect to the overall range of precipitation values. From the analysis of precipitation patterns (Sect. 3.2), we further know that there are strong large-scale precipitation gradients in the regional domain. In the variogram analysis small- and intermediate-scale structures may be hidden by the large-scale gradient. To avoid this and non-stationarity of patterns, we first present variograms of detrended precipitation fields (Sect. 3.4.2). However, to assess processes acting at different scales, variograms of non-detrended precipitation patterns are subsequently analyzed in a scale analysis (Sect. 3.4.3). Finally, a 2-D autocorrelation analysis is used to comment on directional dependencies of precipitation patterns (Sect. 3.4.4).

### 3.4.1 Large-scale precipitation trends

Large-scale precipitation patterns show a strong gradient (Fig. 4). Therefore a plane is fit linearly to the precipitation fields describing the large-scale precipitation trend (Table 2). The trend on 31 January 2016 roughly points toward the north. For the precipitation event on 5 March 2016 the trend points roughly to the south. Given a southerly advection on 5 March 2016 this direction corresponds to the main wind direction and therefore agrees with the theory of large-scale orographic precipitation enhancement or rather the drying trend due to sinking air further downstream within the mountain range. The south–north gradient on 31 January 2016 roughly agrees with the main wind direction but points out that regional trends in larger-scale patterns may not exactly be aligned with wind direction. For the precipitation event on 4 February 2016 the intensity of the trend (i.e., the strength of inclination of the linearly fitted plane) is, however, weak and therefore the orientation of the slope is arbitrary. For this day, we hypothesize that either dynamics were more variable, preventing the evolution of a strong gradient, or lifting condensation due to the orography was not as efficient as for the other two events. For two events (31 January and 4 February 2016) the model has trouble reproducing the trend (i.e., orientation and intensity of the linearly fitted plane). For 31 January 2016 the deviation of orientation between the trends of radar precipitation and WRF precipitation at 2830 m a.s.l. is about 70 but with a similar intensity of the trend. For 4 February 2016 the model shows a strong trend in precipitation, while the intensity of the trend is weak in the entirely filtered radar data. However, for the precipitation event on 5 March 2016 the trend is reasonably captured by the model with a deviation of the orientation of 16.6 and a slightly stronger intensity of the trend in the model compared to the radar estimation.

Disagreement in precipitation patterns, trend orientation and intensity on 4 February 2016 (quite homogeneous precipitation distribution in the radar estimate (Fig. 4) compared to the strong downstream shift of precipitation in WRF) and the overestimation of precipitation in the model give evidence for a too simplistic representation of precipitation in the model (i.e., simplified microphysics and cloud dynamics), which tends to overestimate the effect of the highest topographic features but misses precipitation over shallower areas. Good agreement in the intensity of the trend on 31 January 2016 and good agreement of the orientation of the trend on 5 March 2016, however, show that the model is able to capture large-scale precipitation trends, which may be connected to a large-scale orographic enhancement.

### 3.4.2 Spatial variability in detrended precipitation fields

On the smallest scales a strong difference is visible in variograms of detrended entirely filtered and detrended partly filtered radar precipitation, with weaker variability for entirely filtered data (Fig. 6). The smallest-scale structures in the radar data are likely an indicator of residual noise in the partly filtered radar data (Sect. 2.3). However, it could also imply microscale precipitation features. This stresses the challenge of processing high-resolution radar data (Sect. 2.3) to obtain a reasonable radar precipitation field. In any case the entirely filtered radar precipitation estimates may be regarded as clean concerning residual clutter and will therefore be used for all subsequent analysis.

Figure 6Variograms of detrended snow precipitation normalized by the domain-wide variance of precipitation for the precipitation events on (a) 31 January 2016, (b) 4 February 2016 and (c) 5 March 2016 for the regional domain (450 m horizontal grid spacing, Fig. 1). Variograms are given for partly filtered (red) and entirely filtered (orange) radar snow precipitation, WRF snow precipitation at 2830 m above sea level (m a.s.l., blue) and WRF total ground precipitation (violet). WRF precipitation is from simulations with weak terrain smoothing (Sect. 2.1). All precipitation fields are masked.

Variograms of entirely filtered and detrended radar precipitation show a steep increase in variability on the smallest scales, while the increase in variability becomes weaker for larger scales (less steep slope in the variograms). Small-scale patterns are likely driven by small-scale precipitation cells induced by local cloud dynamics and microphysics. Such small-scale structures are repeated on intermediate scales and lead to a weaker increase in variability, as fewer new spatial features are added. At larger scales variability reaches the total variability in the detrended data.

Compared to radar precipitation WRF precipitation at 2830 m a.s.l. shows a lower variability and a flatter increase in variability at small scales, giving evidence for a smoother precipitation distribution at the smallest scales. The lack of small-scale patterns shows that the radar sees more variability at the smallest scales, while WRF likely misses the smallest-scale processes. Variability in radar and WRF precipitation at 2830 m a.s.l. at large scales (>5 km), especially on 4 February 2016, shows less systematic differences than at small scales. This indicates that, with respect to total variability, patterns at these scales are well represented. Total ground precipitation shows a higher variability compared to precipitation at 2830 m a.s.l. (except for 4 February 2016), which is an indication that near-surface processes are active in the model.

Variograms of precipitation in the local domain (300 m resolution, Fig. 1) look similar to variograms of the regional domain (450 m resolution) but reach domain-wide variability at about 5 km lag distance (Supplement S2), while on the regional scale the domain-wide variability is reached at a distance of about 15–20 km (Fig. 6). Furthermore, the difference between radar and WRF precipitation variability at small scales is larger on the local domain compared to the regional domain. This and a systematic underestimation of precipitation variability at scales <5 km (on the regional domain) compared to the precipitation variability in radar estimates indicate that mountain-ridge-scale precipitation processes are underrepresented in the model.

### 3.4.3 Scale breaks and dominating processes

Scale breaks were previously found to be connected to changes in dominant processes (e.g., Deems et al.2006). Here, we present a scale analysis including variability due to large-scale precipitation processes. Therefore, we present variograms of non-detrended precipitation fields, being aware that a certain portion of the small- and intermediate-scale precipitation variability may become hidden. As precipitation patterns are known to be driven by topography and wind, we present variograms of topography together with the variograms of precipitation.

Variograms of topography clearly reveal two scale breaks (Fig. 7). The first scale break is between 1 and 2.5 km depending on the resolution; the second scale break is at 5 and 6 km for real topography and weakly smoothed WRF topography, respectively. For topography the two scale breaks separate the mountain-slope scale (<1–2 km), mountain-ridge-to-valley scale (between ∼1–2 km and ∼5 km), and the scale of repeated mountain ridges and valleys (>5 km).

Figure 7Normalized variograms of the snow precipitation events on (a) 31 January 2016, (b) 4 February 2016 and (c) 5 March 2016 for the regional domain (450 m horizontal grid spacing, Fig. 1). Variograms are given for entirely filtered radar snow precipitation (orange), WRF snow precipitation at 2830 m above sea level (m a.s.l., blue) and WRF total ground precipitation (violet). Additionally, variograms are given for real topography (based on dhm25 © 2018 swisstopo (5 740 000 000), black) and WRF topography (gray). WRF topography and precipitation are from simulations with weak terrain smoothing (Sect. 2.1). All precipitation fields are masked.

For consistency reasons, all variograms in Fig. 7 are presented with two scale breaks. Scale breaks for all events and both resolutions are basically grouped in two areas (∼1–2.5 and 5–10 km for 450 m resolution, Fig. 7; and ∼800 m–1.2 km and 2.5–5 km for 300 m resolution, Supplement S2), even though for precipitation some scale breaks are arbitrary. Albeit scale breaks of precipitation do not exactly match scale breaks of topography, breaks at similar scales as well as similar slopes of topography and precipitation at small scales support the interpretation of topography-dependent precipitation patterns. On the smallest scales (<1–2 km) the slopes of precipitation variograms are similar to the slopes of the variograms of corresponding topography. This is an indication that precipitation patterns on mountain-slope scales may be terrain driven. Processes acting at these scales could be small-scale cloud dynamical processes such as the seeder–feeder mechanism or preferential deposition . The latter is, however, for most mountain ridges unlikely to be seen in precipitation fields at 2830 m a.s.l. as it happens close to the ground. For the precipitation event on 5 March 2016, on scales >5–7 km (i.e., for the scales above the second scale break), the slopes of the normalized variograms of radar and WRF precipitation at radar elevation are similar. Large-scale gradients at these scales are most likely driven by large-scale orographic precipitation enhancement (e.g., Stoelinga et al.2013). Good agreement of the slopes in normalized variograms between radar and WRF precipitation is an indicator that the model has the potential to properly represent the strength of the large-scale gradient with respect to the overall variability, i.e., large-scale orographic precipitation enhancement. Disagreement of variograms of precipitation and topography at these scales further supports the hypothesis that largest-scale precipitation is mainly determined by large-scale orographic precipitation enhancement, which introduces an increase in variability in precipitation at large scales, while large-scale topography reveals a repeated pattern of valleys and peaks (i.e., constant variability). Overall, this analysis supports the hypothesis in Sect. 3.2 that precipitation patterns in the regional domain are topography driven.

### 3.4.4 Two-dimensional variability patterns

Finally, the combined influence of topography and the general wind direction on snow precipitation patterns in the regional domain is assessed by spatial 2-D autocorrelation maps (Fig. 8). Like variograms, autocorrelation is dependent on large-scale trends. The general direction of 2-D autocorrelation patterns is the same for detrended (Fig. 8) and non-detrended (not shown) precipitation patterns. However, autocorrelation patterns of detrended precipitation fields show much shorter decorrelation lengths. This is due to the spatial coherence introduced by large-scale trends in precipitation. To avoid biased autocorrelation data, only 2-D autocorrelation maps of detrended precipitation fields are shown. However, we keep in mind that large-scale trends are present.

Autocorrelation maps of topography (Fig. 8a and e) represent a northwest-to-southeast-oriented pattern, which is, although weaker, repeated in the southwest-to-northeast and west-to-east directions. For snow precipitation events with dominating northwesterly to northerly advection, the main axis of the snow precipitation 2-D autocorrelation pattern is oriented in a northwest-to-southeast direction and therefore in alignment with both topography and the main wind direction (Fig. 8b–c and f–g). Patterns of WRF precipitation at 2830 m a.s.l. are rotated toward a north–south direction on 4 February 2016. For dominating southerly advection the 2-D autocorrelation map of radar precipitation shows a more homogeneous pattern compared to autocorrelation patterns for northern to northwestern advection but a weak southwest-to-northeast orientation of larger-scale patterns (Fig. 8d). For the WRF simulations a strong southwest-to-northeast orientation is present in the autocorrelation map for the precipitation event on 5 March 2016 (Fig. 8h). Even though isotropic variograms reveal good agreement in domain-wide variability, 2-D autocorrelation maps show that this may not necessarily go along with good agreement of the orientation of patterns. The best agreement in the orientation of patterns is found for 31 January 2016. For the three events, 2-D autocorrelation maps of detrended precipitation reveal a smoother distribution of precipitation on the smallest scales in the model compared to radar data due to fewer small-scale structures in the model. Conversely, a strong decrease in autocorrelation in the east–west direction is visible for 5 March 2016. This shows that WRF simulations have a stronger dependency on both wind direction and topography and tend to generate strong precipitation bands in the main wind direction, confirming the overly simplistic behavior of the model.

For ground precipitation, 2-D autocorrelation patterns tend to be repeated in the southwest-to-northeast and west-to-east directions as seen for topography (Fig. 8j–g). This stresses the hypothesis that the influence of topographic features on WRF ground precipitation is stronger than at radar elevation and gives evidence that these results are likely produced by near-surface topographically driven pre-depositional processes such as preferential deposition or the seeder–feeder mechanism in the model, for example. While a topography dependency was already found in isotropic variograms, this 2-D autocorrelation analysis reveals that the wind direction additionally strongly impacts the snow precipitation distribution.

Figure 8Spatial 2-D autocorrelation maps for the regional domain (450 m resolution) of detrended (a) real topography (based on dhm25 © 2018 swisstopo (5740 000 000)), (b–d) entirely filtered radar snow precipitation, (e) WRF topography, (f–h) WRF snow precipitation at 2830 m above sea level (m a.s.l.) and (j–l) WRF total ground precipitation. Autocorrelation maps of snow precipitation are for the three snow precipitation events on 31 January, 4 February and 5 March 2016. WRF topography and precipitation are from simulations with weak terrain smoothing (Sect. 2.1). Radar precipitation and WRF precipitation at 2830 m a.s.l. are masked (as shown in Fig. 4).

## 3.5 Dependence of spatial variability on model resolution and smoothing

Geostatistical analyses presented in this study demonstrate that precipitation on the regional scale (>5 km) is reasonably represented in the WRF model, while small-scale precipitation variability is systematically underestimated in the model simulations with a horizontal grid spacing of 450 m (Sect. 3.4). Variograms up to a maximum lag distance of 5 km for domain Dischma (Fig. 1) reveal an increase in variability for increasing model resolution (Fig. 9). However, simulated variability stays far below the variability in entirely filtered radar precipitation. Depending on the event an increase in variability is present for 150 and 50 m resolution. This indicates that the smallest-scale precipitation dynamics are still not fully resolved at 50 m resolution. A comparison of variograms for simulations with strongly smoothed terrain compared to simulations with weaker terrain smoothing (Sect. 2.1) reveal that a stronger terrain smoothing may result in less explained variability in normalized variograms (not shown). Even though this signal is not consistent for all events, we can show that a better representation of topography due to higher resolution and less smoothing has the potential to increase the explained variability in precipitation patterns. An increase in variability at small scales (<5 km) indicates that more small-scale patterns are resolved at higher resolutions (50 m horizontal grid spacing) in the model. Our simulations are currently limited to the presented resolutions and strong terrain smoothing due to model instabilities. However, based on the presented results, once available, the immersed boundary method version of WRF will likely be a good tool to allow for steeper slopes in the simulation and moving toward higher-resolution LES simulations to resolve further small-scale wind fields, which drive the precipitation structures.

Figure 9Variograms of detrended snow precipitation normalized by the domain-wide variance of precipitation for the precipitation events on (a) 31 January 2016, (b) 4 February 2016 and (c) 5 March 2016 for domain Dischma (Fig. 1). Variograms are given for entirely filtered radar snow precipitation (red) and WRF snow precipitation at 2830 m above sea level (m a.s.l.) with 450 m (violet), 150 m (blue) and 50 m (light blue) resolution. WRF snow precipitation is from simulations with weak terrain smoothing (Sect. 2.1). Radar precipitation is masked.

4 Conclusions and outlook

The implementation of COSMO–WRF is a further step in performing very-high-resolution precipitation simulations in complex alpine terrain to address the question of the relative importance of cloud dynamics and particle–flow interactions on a mountain-ridge scale. In this validation study, we show that COSMO–WRF is able to reasonably simulate temperature, relative humidity and wind direction but tends to (strongly) overestimate near-surface wind speeds. This may be due to many reasons from an overestimation of speedup effects on the windward side of mountain ridges to an underrepresentation of small terrain features. Additional reasons are likely but remain unknown and future work is needed to address these issues. Relative humidity patterns are highly variable and may be a sign that subsidence and lifting produce too strong effects in the (partially parameterized) cloud dynamics, given the good representation of topography at larger scales.

Regional- and local-scale precipitation patterns in the COSMO–WRF simulations are in partially good agreement with MeteoSwiss operational radar measurements and automatic weather stations. For the three events analyzed here, precipitation estimates from WRF simulations are higher compared to precipitation estimates from radar measurements. A general overestimation of precipitation produced by WRF is consistent with an overestimation of subsidence and lifting. Overestimation of precipitation in WRF simulations has been documented previously for snow precipitation over complex terrain (e.g., Silverman et al.2013); among other reasons it may be due to the high model resolution and therefore more complex topography and higher mountain peaks compared to common high-resolution simulations.

To specifically address processes such as the seeder–feeder mechanism or preferential deposition, an analysis of hydrometeors and precipitation distributions in vertical profiles across mountain ridges is needed. To connect pre-depositional processes with post-depositional processes, even higher-resolution WRF simulations would be required. This might be achieved by employing the immersed-boundary method version setup of WRF. A parameterization of post-depositional processes in WRF or using WRF simulations as a boundary condition for simulations with the Alpine surface processes model Alpine3D would then allow validation of modeled snow accumulation patterns compared to measured snow accumulation patterns. Furthermore, simulations of precipitation patterns in complex terrain need to be analyzed with higher temporal resolution (e.g., on the order of minutes), as contributing processes show high temporal variability. Future work will include addressing the temporal variability in precipitation patterns using radar observations, along with an analysis of precipitation growth with respect to topography and wind direction.

Code and data availability
Code and data availability.

A COSMO–WRF documentation is published in . Radar and COSMO data can be made available upon request.

Appendix A: Morrison microphysics in WRF

The Morrison microphysics scheme includes prognostic equations of number concentration and mass mixing ratio of five precipitation species (rain, snow, ice, graupel and cloud droplets). The parametrization of rain, snow, ice and cloud droplets is based on . The implementation of graupel follows , except for minimum mixing ratios, which are required to produce graupel from the collision of rain and snow, snow and cloud water, and rain and cloud ice, which are based on .

The kinetic equations include advection, sedimentation and turbulent diffusion as well as source and sink terms of ice nucleation and droplet activation, condensation and deposition, coalescence and diffusional growth, collection, melting and freezing, and ice multiplication . For graupel deposition, collection, collision, accretion, freezing and melting processes are parameterized .

Size distribution functions are gamma functions:

$\begin{array}{}\text{(A1)}& N\left(D\right)={N}_{\mathrm{0}}{D}^{\mathit{\mu }}{e}^{-\mathit{\lambda }D},\end{array}$

where D is the particle diameter and μ is the shape parameter of the distribution function, which is μ=0 for rain, snow, ice and graupel, resulting in an exponential function for N(D). λ and N0 are the slope and intercept, respectively, of the size distribution, evaluated by the predicted number concentration N and mass mixing ratio q:

$\begin{array}{}\text{(A2)}& \mathit{\lambda }={\left[\frac{cN\mathrm{\Gamma }\left(\mathit{\mu }+d+\mathrm{1}\right)}{q\mathrm{\Gamma }\left(\mathit{\mu }+\mathrm{1}\right)}\right]}^{\mathrm{1}/d}\end{array}$

and

$\begin{array}{}\text{(A3)}& {N}_{\mathrm{0}}=\frac{N{\mathit{\lambda }}^{\mathit{\mu }+\mathrm{1}}}{\mathrm{\Gamma }\left(\mathit{\mu }+\mathrm{1}\right)},\end{array}$

where Γ is the gamma function. c and d are the parameters of the power-law function m=cDd indicating the mass–diameter relationship. Terminal fall speeds are also assumed to have a power-law form of $v\left(D\right)=\frac{{\mathit{\rho }}_{\text{sur}}}{\mathit{\rho }}a{D}^{b}$, with individual parameters a and b for the different species. ρ is the air density and ρsur the air density at sea level. For simplification all species are assumed to be spheres. Additionally, the particles do not have any particle inertia.

Author contributions
Author contributions.

FG and NB performed the analysis, which was supported by fruitful discussions with ML, RM, AB and UG. FG, VS, MD, RM and ML developed the COSMO–WRF coupling and simulation setups, which were run by FG. NB and MG processed the radar data. FG, ML, RM, NB and AB contributed to the design of the concept. FG and NB, with contributions from all authors, prepared the paper.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

Edited by: Ross Brown
Reviewed by: Heini Wernli and two anonymous referees

References

Arnold, D., Schicker, I., and Seibert, P.: High-Resolution Atmospheric Modelling in Complex Terrain for Future Climate Simulations(HiRmod), Report 2010, Tech. rep., Institute of Meteorology (BOKU-Met), University of Natural Resources and Life Sciences, Vienna, Austria, 2010. a

Arthur, R., Lundquist, K. A., Mirocha, J. D., Hoch, S. W., and Chow, F. K.: High-resolution simulations of downslope flows over complex terrain using WRF-IBM, 17th Conference on Mountain Meteorology, American Meteorological Society, Paper 7.6, 18 pp., 2016. a

Beljaars, A. C. M.: The parameterization of surface fluxes in large-scale models under free convection, Q. J. Roy. Meteor. Soc., 121, 255–270, https://doi.org/10.1002/qj.49712152203, 1994. a

Bergeron, T.: On the low-level redistribution of atmospheric water caused by orography, Suppl. Proc. Int. Conf. Cloud Phys., Tokyo, 96–100, 1965. a, b

Besic, N., Figueras i Ventura, J., Grazioli, J., Gabella, M., Germann, U., and Berne, A.: Hydrometeor classification through statistical clustering of polarimetric radar measurements: a semi-supervised approach, Atmos. Meas. Tech., 9, 4425–4445, https://doi.org/10.5194/amt-9-4425-2016, 2016. a

Caldwell, P., Chin, H., Bader, D., and Bala, G.: Evaluation of a WRF dynamical downscaling simulation over California, Clim. Change, 95, 499–521, https://doi.org/10.1007/s10584-009-9583-5, 2009. a

Choularton, T. W. and Perry, S. J.: A model of the orographic enhancement of snowfall by the seeder-feeder mechanism, Q. J. Roy. Meteor. Soc., 112, 335–345, https://doi.org/10.1002/qj.49711247204, 1986. a

Colle, B.: Sensitivity of orographic precipitation to changing ambient conditions and terrain geometries: An idealized modelling perspective, J. Atmos. Sci., 61, 588–606, https://doi.org/10.1175/1520-0469(2004)061<0588:SOOPTC>2.0CO;2, 2004. a

Cooper, S. J., Wood, N. B., and L'Ecuyer, T. S.: A variational technique to estimate snowfall rate from coincident radar, snowflake, and fall-speed observations, Atmos. Meas. Tech., 10, 2557–2571, https://doi.org/10.5194/amt-10-2557-2017, 2017. a

Dadic, R., Mott, R., Lehning, M., and Burlando, P.: Wind influence on snow depth distribution and accumulation over glaciers, J. Geophys. Res., 115, F01012, https://doi.org/10.1029/2009JF001261, 2010. a

Daniels, M. H., Lundquist, K. A., Mirocha, J. D., Wiersema, D. J., and Chow, F. K.: A New Vertical Grid Nesting Capability in the Weather Research and Forecasting (WRF) Model, Mon. Weather Rev., 144, 3725–3747, https://doi.org/10.1175/MWR-D-16-0049.1, 2016. a

Deems, J. S., Fassnacht, S. R., and Elder, K. J.: Fractal Distribution of Snow Depth from Lidar Data, J. Hydrometeor., 7, 285–297, https://doi.org/10.1175/JHM487.1, 2006. a, b, c, d

Deems, J. S., Fassnacht, S. R., and Elder, K. J.: Interannual Consistency in Fractal Snow Depth Patterns at Two Colorado Mountain Sites, J. Hydrometeor., 9, 977–988, https://doi.org/10.1175/2008JHM901.1, 2008. a, b

Dore, A. J., Choularton, T. W., Fowler, D., and Crossely, A.: Orographic enhancement of snowfall, Environ. Pollut., 75, 175–179, https://doi.org/10.1016/0269-7491(92)90037-B, 1992. a

Dyer, A. J. and Hicks, B. B.: Flux-gradient relationships in the constant flux layer, Q. J. Roy. Meteor. Soc., 96, 715–721, https://doi.org/10.1002/qj.49709641012, 1970. a

European Environmental Agency: CORINE Land Cover (CLC) 2006 raster data, Version 13, 2006. a

Gabella, M., Speirs, P., Hamann, U., Germann, U., and Berne, A.: Measurement of Precipitation in the Alps Using Dual-Polarization C-Band Ground-Based Radars, the GPM Spaceborne Ku-Band Radar, and Rain Gauges, Remote Sens., 9, 1147, https://doi.org/10.3390/rs9111147, 2017. a

Gerber, F. and Sharma, V.: Running COMO-WRF on very high resolution over complex terrain, Laboratory of Cryospheric Sciences, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, https://doi.org/10.16904/envidat.35, 2018. a, b, c, d

Gerber, F., Lehning, M., Hoch, S. W., and Mott, R.: A close-ridge small-scale atmospheric flow field and its influence on snow accumulation, J. Geophys. Res.-Atmos., 122, 7737–7754, https://doi.org/10.1002/2016JD026258, 2017. a, b

Gerber, F., Sharma, V., Mott, R., Daniels, M. and Lehning, M.: DISCHMEX – High-resolution WRF simulations in complex alpine terrain and station measurements, Laboratory of Cryospheric Sciences, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, https://doi.org/https://doi.org/10.16904/envidat.50, 2018.

Germann, U., Galli, G., Boscacci, M., and Bolliger, M.: Radar precipitation measurement in a mountainous region, Q. J. Roy. Meteor. Soc., 132, 1669–1692, https://doi.org/10.1256/qj.05.190, 2006. a

Germann, U., Boscacci, M., Gabella, M., and Sartori, M.: Peak performance: Radar design for prediction in the Swiss Alps, Meteorological Technology International, 42–45, 2015. a, b

Goger, B., Rotach, M. W., Gohm, A., Fuhrer, O., Stiperski, I., and Holtslag, A. M.: The Impact of Three-Dimensional Effects on the Simulation of Turbulence Kinetic Energy in a Major Alpine Valley, Bound.-Lay. Meteorl., 168, 1–27, https://doi.org/10.1007/s10546-018-0341-y, 2018. a

Gómez-Navarro, J. J., Raible, C. C., and Dierer, S.: Sensitivity of the WRF model to PBL parametrisations and nesting techniques: evaluation of wind storms over complex terrain, Geosci. Model Dev., 8, 3349–3363, https://doi.org/10.5194/gmd-8-3349-2015, 2015. a, b, c

Grünewald, T., Schirmer, M., Mott, R., and Lehning, M.: Spatial and temporal variability of snow depth and ablation rates in a small mountain catchment, The Cryosphere, 4, 215–225, https://doi.org/10.5194/tc-4-215-2010, 2010. a

Grünewald, T., S., D., Cattin, R., Steiner, P., Steinkogler, W., Fundel, F., and Lehning, M.: Mapping frequencies of icing on structures in Switzerland, J. Eng. Ind. Aerodynam., 107–108, 76–82, https://doi.org/10.1016/j.jweia.2012.03.022, 2012. a, b

Hong, S.-Y., Noh, Y., and Dudhia, J.: A New Diffusion Package with an Explicit Treatment of Entrainment Processes, Mon. Weather Rev., 134, 2318–2341, https://doi.org/10.1175/MWR3199.1, 2006. a

Houze Jr., R. A.: Orographic effects on precipitating clouds, Rev. Geophys., 50, RG1001, https://doi.org/10.1029/2011RG000365, 2012. a, b

Huwald, H., Higgins, C. W., Boldi, M.-O., Bou-Zeid, E., Lehning, M., and Parlange, M. B.: Albedo effect on radiative errors in air temperature measurements, Water Resour. Res., 45, W08431, https://doi.org/10.1029/2008WR007600, 2009. a, b

Jiménez, P. A. and Dudhia, J.: Improving the Representation of Resolved and Unresolved Topographic Effects on Surface Wind in the WRF Model, J. Appl. Meteorol. Climatol., 51, 300–316, https://doi.org/10.1175/JAMC-D-11-084.1, 2012. a, b, c

Lehning, M. and Fierz, C.: Assessment of snow transport in avalanche terrain, Cold Reg. Sci. Technol., 51, 240–252, https://doi.org/10.1016/j.coldregions.2007.05.012, 2008. a

Lehning, M., Löwe, H., Ryser, M., and Raderschall, N.: Inhomogeneous precipitation distribution and snow transport in steep terrain, Water Resour. Res., 44, W07404, https://doi.org/10.1029/2007WR006545, 2008. a, b, c, d

Leung, L. R. and Qian, Y.: The sensitivity of precipitation and snowpack simulations to model resolution via nesting in regions of complex terrain, J. Hydrometeor., 4, 1025–1043, https://doi.org/10.1175/1525-7541(2003)004<1025:TSOPAS>2.0.CO;2, 2003. a, b

Liu, C., Ikeda, K., Thompson, G., Rasmussen, R., and Dudhia, J.: High-Resolution Simulations of Wintertime Precipitation in the Colorado Headwaters Region: Sensitivity to Physics Parameterizations, Mon. Weather Rev., 139, 3533–3553, https://doi.org/10.1175/MWR-D-11-00009.1, 2011. a

Lundquist, K. A., Chow, F. K., and Lundquist, J. K.: An immersed boundary method for the Weather Research and Forecasting model, Mon. Weather Rev., 138, 796–817, https://doi.org/10.1175/2009MWR2990.1, 2010. a

Lundquist, K. A., Chow, F. K., and Lundquist, J. K.: An Immersed Boundary Method Enabling Large-Eddy Simulations of Flow over Complex Terrain in the WRF Model, Mon. Weather Rev., 140, 3936–3955, https://doi.org/10.1175/MWR-D-11-00311.1, 2012. a

Ma, Y. and Liu, H.: Large-Eddy Simulations of Atmospheric Flows Over Complex Terrain Using the Immersed-Boundary Method in the Weather Research and Forecasting Model, Bound.-Lay. Meteorol., 165, 421–445, https://doi.org/10.1007/s10546-017-0283-9, 2017. a

Mass, C., Ovens, D., Westrick, K., and Colle, B. A.: Does increasing horizontal resolution produce more skillful forecasts?, B. Am. Meteorol. Soc., 83, 407–430, https://doi.org/10.1175/1520-0477(2002)083<0407:DIHRPM>2.3.CO;2, 2002. a, b

METI/NASA: 2009, ASTER Global Digital Elevation Model V002, NASA EOSDIS Land Processes DAAC, USGS Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, https://doi.org/10.5067/ASTER/ASTGTM.002, 2009. a

Mirocha, J., Kosović, B., and Kirkil, G.: Resolved Turbulence Characteristics in Large-Eddy Simulations Nested within Mesoscale Simulations Using the Weather Research and Forecasting Model, Mon. Weather Rev., 142, 806–831, https://doi.org/10.1175/MWR-D-13-00064.1, 2014. a

Mirocha, J. D., Lundquist, J. K., and Kosović, B.: Implementation of a Nonlinear Subfilter Turbulence Stress Model for Large-Eddy Simulation in the Advanced Research WRF Model, Mon. Weather Rev., 138, 4212–4228, https://doi.org/10.1175/2010MWR3286.1, 2010. a

Morrison, H., Curry, J. A., and Khvorostyanov, V. I.: A New Double-Moment Microphysics Parameterization for Application in Cloud and Climate Models, Part I: Description, J. Atmos. Sci., 62, 1665–1677, https://doi.org/10.1175/JAS3446.1, 2005. a, b, c, d

Morrison, H., Thompson, G., and Tatarskii, V.: Impact of Cloud Microphysics on the Development of Trailing Stratiform Precipitation in a Simulated Squall Line: Comparison of One- and Two-Moment Schemes, Mon. Weather Rev., 137, 991–1007, https://doi.org/10.1175/2008MWR2556.1, 2009. a

Mott, R. and Lehning, M.: Meteorological Modeling of Very High-Resolution Wind Fields and Snow Deposition for Mountains, J. Hydrometeor., 11, 934–949, https://doi.org/10.1175/2010JHM1216.1, 2010. a

Mott, R., Schirmer, M., Bavay, M., Grünewald, T., and Lehning, M.: Understanding snow-transport processes shaping the mountain snow-cover, The Cryosphere, 4, 545–559, https://doi.org/10.5194/tc-4-545-2010, 2010. a, b, c, d

Mott, R., Schirmer, M., and Lehning, M.: Scaling properties of wind and snow depth distribution in an Alpine catchment, J. Geophys. Res., 116, D06106, https://doi.org/10.1029/2010JD014886, 2011. a, b

Mott, R., Scipión, D., Schneebeli, M., Dawes, N., Berne, A., and Lehning, M.: Orographic effects on snow deposition patterns in mountainous terrain, J. Geophys. Res.-Atmos., 119, 1419–1439, https://doi.org/10.1002/2013JD019880, 2014. a, b, c, d

Mott, R., Daniels, M., and Lehning, M.: Atmospheric Flow Development and Associated Changes in Turbulent Sensible Heat Flux over Patchy Mountain Snow Cover, J. Hydrometeor., 16, 1315–1340, https://doi.org/10.1175/JHM-D-14-0036.1, 2015. a

Mott, R., Schlögl, S., Dirks, L., and Lehning, M.: Impact of Extreme Land Surface Heterogeneity on Micrometeorology over Spring Snow Cover, J. Hydrometeor., 18, 2705–2722, https://doi.org/10.1175/JHM-D-17-0074.1, 2017. a

Muñoz Esparza, D., Lundquist, J. K., Sauer, J. A., Kosovic, B., and Linn, R. R.: Coupled mesoscale-LES modeling of a diurnal cycle during the CWEX-13 field campaign: From weather to boundary-layer eddies, J. Adv. Model. Earth Sy., 9, 1572–1594, https://doi.org/10.1002/2017MS000960, 2017. a

Niu, G.-Y., Yang, Z.-L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Xia, Y.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res., 116, D12109, https://doi.org/10.1029/2010JD015139, 2011. a

Obukhov, A. M.: Turbulence in an atmosphere with a non-uniform temperature, Bound.-Lay. Meteorol., 2, 7–29, 1971. a

Panziera, L., James, C. N., and Germann, U.: Mesoscale organization and structure of orographic precipitation producing flash floods in the Lago Maggiore region, Q. J. Roy. Meteor. Soc., 141, 224–248, https://doi.org/10.1002/qj.2351, 2015. a

Paulson, C. A.: The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer, J. Appl. Meteorol., 9, 857–861, https://doi.org/10.1175/1520-0450(1970)009<0857:TMROWS>2.0.CO;2, 1970. a

Pineda, N., Jorba, O., Jorge, J., and Baldasano, J. M.: Using NOAA AVHRR and SPOT VGT data to estimate surface parameters: application to a mesoscale meteorological model, Int. J. Remote Sens., 25, 129–143, https://doi.org/10.1080/0143116031000115201, 2004. a

Pohl, B.: Testing WRF capability in simulating the atmospheric water cycle over Equatorial East Africa, Clim. Dynam., 37, 1375–1379, https://doi.org/10.1007/s00382-011-1024-2, 2011. a

Prokop, A.: Assessing the applicability of terrestrial laser scanning for spatial snow depth measurements, Cold Reg. Sci. Technol., 54, 155–163, https://doi.org/10.1016/j.coldregions.2008.07.002, 2008. a

Purdy, J. C., Austin, G. L., W., S. A., and Cluckie, I. D.: Radar evidence of orographic enhancement due to the seeder feeder mechanism, Meteorol. Appl., 12, 199–206, https://doi.org/10.1017/S1350482705001672, 2005. a, b

Reisner, J., Rasmussen, R. M., and Bruintjes, R. T.: Explicit forecasting of supercooled liquid water in winter storms using the MM5 mesoscale model, Q. J. Roy. Meteor. Soc., 124, 1071–1107, https://doi.org/10.1002/qj.49712454804, 1998. a, b

Rutledge, S. A. and Hobbs, P. V.: The mesoscale and microscale structure of organization of clouds and precipitation in midlatitude cyclones. XII: A diagnostic modeling study of precipitation development in narrow cold-frontal rainbands, J. Atmos. Sci., 41, 2949–2972, https://doi.org/10.1175/1520-0469(1984)041<2949:TMAMSA>2.0.CO;2, 1984. a

Saltikoff, E., Lopez, P., Taskinen, A., and Pulkkinen, S.: Comparison of quantitative snowfall estimates from weather radar, rain gauges and a numerical weather prediction model, Boreal Environ. Res., 20, 667–678, 2015. a

Schirmer, M. and Lehning, M.: Persistence in intra-annual snow depth distribution: 2. Fractal analysis of snow depth development, Water Resour. Res., 47, W09517, https://doi.org/10.1029/2010WR009429, 2011. a, b

Schirmer, M., Wirz, V., Clifton, A., and Lehning, M.: Persistence in intra-annual snow depth distribution: 1. Measurements and topographic control, Water Resour. Res., 47, W09516, https://doi.org/10.1029/2010WR009426, 2011. a, b, c, d

Schlögl, S., Lehning, M., and Mott, R.: Representation of horizontal transport processes in snowmelt modelling by applying a footprint approach, Front. Earth Sci., https://doi.org/10.3389/feart.2018.00120, online first, 2018. a

Schmucki, E., Marty, C., Fierz, C., Weingartner, R., and Lehning, M.: Impact of climate change in Switzerland on socioeconomic snow indices, Theor. Appl. Climatol., 127, 875–889, https://doi.org/10.1007/s00704-015-1676-7, 2017. a

Scipión, D. E., Mott, R., Lehning, M., Schneebeli, M., and Berne, A.: Seasonal small-scale spatial variability in alpine snowfall and snow accumulation, Water Resour. Res., 49, 1446–1457, https://doi.org/10.1002/wrcr.20135, 2013. a, b, c

Silverman, N. L., Maneta, M. P., Chen, S.-H., and Harper, J. T.: Dynamically downscaled winter precipitation over complex terrain of the Central Rockies of Western Montana, USA, Water Resour. Res., 49, 458–470, https://doi.org/10.1029/2012WR012874, 2013. a, b, c, d, e

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Duda, M. G., Huang, X.-Y., Wang, W., and Powers, J. G.: A Description of the Advanced Research WRF Version 3, Tech. rep., Mesoscale and Microscale Meteorological Division, National Center for Atmospheric Research, Boulder, Colorado, USA, 2008. a, b

Stoelinga, M. T., Stewart, R. E., Thompson, G., and Thériault, J. M.: Mountain Weather Research and Forecasting: Recent Progress and Current Challenges, in: Microphysical Processes Within Winter Orographic Cloud and Precipitation Systems, Springer Netherlands, 345–408, https://doi.org/10.1007/978-94-007-4098-3_7, 2013. a, b, c, d

Talbot, C., Bou-Zeid, E., and Smith, J.: Nested Mesoscale Large-Eddy Simulations with WRF: Performance in Real Test Cases, J. Hydrometeor., 13, 1421–1441, https://doi.org/10.1175/JHM-D-11-048.1, 2012. a

Tedesche, M. E., Fassnacht, S. R., and Meiman, P. J.: Scales of snow depth variability in high elevation rangeland sagebush, Front. Earth Sci., 11, 469–481, https://doi.org/10.1007/s11707-017-0662-z, 2017. a

Trujillo, E., Molotch, N. P., Goulden, M. L., Kelly, A. E., and Bales, R. C.: Elevation-dependent influence of snow accumulation on forest greening, Nat. Geosci., 5, 705–709, https://doi.org/10.1038/ngeo1571, 2012. a

Vionnet, V., Martin, E., Masson, V., Lac, C., Naaim Bouvet, F., and Guyomarc'h, G.: High-resolution large eddy simulation of snow accumulation in alpine terrain, J. Geophys. Res.–Atmos., 122, 11005–11021, https://doi.org/10.1002/2017JD026947, 2017. a, b, c, d

Webb, E. K.: Profile relationships: The log-linear range, and extension to strong stability, Q. J. Roy. Meteor. Soc., 96, 67–90, https://doi.org/10.1002/qj.49709640708, 1970. a

Wyngaard, J. C.: Toward Numerical Modeling in the “Terra Incognita”, J. Atmos. Sci., 61, 1816–1826, https://doi.org/10.1175/1520-0469(2004)061<1816:TNMITT>2.0.CO;2, 2004. a

Yang, Z.-L., Niu, G.-Y., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Longuevergne, L., Manning, K., Niyogi, D., Tewari, M., and Xia, Y.: The community Noah land surface model with multiparameterization options (Noah-MP): 2. Evaluation over global river basins, J. Geophys. Res., 116, D12110, https://doi.org/10.1029/2010JD015140, 2011.  a

Zängl, G.: The temperature dependence of small-scale orographic precipitation enhancement, Q. J. Roy. Meteor. Soc., 134, 1167–1181, https://doi.org/10.1002/qj.267, 2008. a, b

Zängl, G., Aulehner, D., Wastl, C., and Pfeiffer, A.: Small-scale precipitation variability in the Alps: Climatology in comparison with semi-idealized numerical simulations, Q. J. Roy. Meteor. Soc., 134, 1865–1880, https://doi.org/10.1002/qj.311, 2008. a

Zhang, D.-L. and Anthes, R. A.: A High-Resolution Model of the Planetary Boundary Layer – Sensitivity Tests and Comparisons with SESAME-79 Data, J. Appl. Meteorol., 21, 1594–1609, https://doi.org/10.1175/1520-0450(1982)021<1594:AHRMOT>2.0.CO;2, 1982. a