Fine-Resolution WRF Simulation of Stably Stratified Flows in Shallow Pre-Alpine Valleys: A Case Study of the KASCADE-2017 Campaign

Abstract

In an overall approach aiming at the development and qualification of various tools designed to diagnose and/or forecast the flows at the local scale in complex terrain, we qualified a numerical model based on the WRF platform and operated in a two-way nested domain mode, down to a horizontal resolution of 111 m for the smallest domain. The area in question is the Cadarache valley (CV), in southeast France, which is surrounded by hills and valleys of various sizes. The CV dimensions (1 km wide and 100 m deep) favor the development of local flows greatly influenced by the diurnal cycle and are prone to thermal stratification, especially during stable conditions. This cycle was well documented due to permanent observations and dedicated field campaigns. These observations were used to evaluate the performance of the model on a specific day among the intensive observation periods carried out during the KASCADE-2017 campaign. The model reproduced the wind flow and its diurnal cycle well, notably at the local CV scale, which constitutes considerable progress with respect to the performances of previous WRF simulations conducted in this area with kilometric resolution, be it operational weather forecasts or dedicated studies conducted on specific days. The diurnal temperature range is underestimated however, together with the stratification intensity of the cold pool observed at night. Consequently, the slope drainage flows along the CV sidewalls are higher in the simulation than in the observations, and the resulting scalar fields (such as specific humidity) are less heterogeneous in the model than in the observations.

1. Introduction

The flows over complex terrain represent a challenge for both their observation and numerical simulation. At a large scale, the major mountain chains constrain the flows in the entire troposphere, and their effect has to be adequately taken into account for an accurate weather forecast. To date, one may consider this question as having been satisfactorily addressed, although some improvements are still required [1]. At the local scale, the topography has a direct impact on the lower layers of the atmosphere (the atmospheric boundary layer; ABL), with a marked diurnal cycle. While at a large scale, the mechanical effects of the topography are dominant, at the local scale, the flow is driven by a combination of both mechanical and buoyancy effects. Furthermore, what we call “local” in fact covers quite a large range of scales, ranging from the smaller hill–valley systems (of O(10 m) in height and O(100 m) in length) to systems of one or two orders of magnitude higher. All these systems interact with each other, but on different time scales.

Such complexity results in great difficulty in monitoring the local transport, and stagnation/blocking effects can generate air quality concerns with polluted air masses trapped in valleys even in rural areas. In cold pool systems, the stagnation can be extreme with the consequence of very low temperatures sometimes observed on the ground (Lehner et al. [2] mentioned a recorded minimum temperature of −52.6 °C observed in the Grünloch basin in 1932). In the French Alps, the Arve valley has recently received a great deal of attention due to dedicated observations and high-resolution numerical simulations [3,4,5]. For a localized potential pollution source, a question arises regarding the behavior of the particle plume, that is, how it is transported and dispersed in a complex terrain, especially during stable stratification conditions with suppressed turbulent vertical mixing. In complex terrain, simple dispersion models fail in the prediction of plume behavior, and the interaction between orographic systems of different forms and sizes may produce unexpected stagnation effects [6,7].

According to the academic interest in local process knowledge, and to the societal importance of the possible impacts of complex topography, many studies have been launched over recent decades, often based on observation field campaigns. We cannot report an exhaustive list of these programs here. We will therefore only mention a few relevant campaigns. Density-driven valley and slope flows and/or cold air pooling were studied during VTMX [8], COLPEX [9], METCRAX [10], MATERHORN [11], ASCOT [12], and the Passy project [13,14]. Rotach et al. [15] published a description of the i-Box setup in the Inn valley. Recently, the Perigão 2017 campaign [16] observed and quite successfully simulated the local valley flows and their reversal during the diurnal cycle. Other studies focused on the effect of slope orientation on surface fluxes and local flows [17,18,19]. All these campaigns saw the deployment of meteorological observation means at the surface and aloft, including remote sensing platforms (sodars, radars, and lidars), among which the scanning Doppler lidar offers a high potential for the fine study of flows in stable layers in which oscillations are often observed (e.g., [20]).

The Cadarache CEA site (CEA stands for “Commissariat à l’Énergie Atomique et aux Énergies Alternatives”, the French nuclear research agency) is located in southeastern France, on the border of the Alps mountain chain and the Rhône river. The area has a complex system of hills and valleys of different sizes (a description is provided in Section 2). The chronic or accidental release of contaminants requires the ability to monitor the transport and dispersion of effluents in this area, especially during stable stratification periods when local, thermally driven valley winds and/or cold pool systems may be established. For a decade, therefore, an effort has been made to better understand the meteorology at a very local scale, and to develop various simulation tools. In particular, the KASCADE-2013 campaign [21] observed the diurnal cycle of valley winds and their stacking in the area, whereas the KASCADE-2017 campaign [22], the observations of which are the base of the present study, focused on the spatial variability of the flow and the stratification inside the valley. Various ways were thus explored to nowcast or forecast the local scale flows. Duine et al. [23] developed a simple but robust method with which to diagnose the presence of a down valley wind in stable conditions based on the observation of a stratification index. Dupuy et al. [24] used artificial neural networks (ANNs) in order to estimate the wind speed and direction in a small valley from routine observations representative of a larger scale. Finally, [25] extended this ANN technique to forecast the valley wind from operational weather simulations performed with a resolution which was too coarse to represent the fine-scale flows.

In this paper, we present an attempt to simulate the flows by running a weather model with a fine (~100 m) horizontal resolution. Recent studies have demonstrated the potential of such simulations in complex areas (e.g., [4,5] for the Arve valley with the Meso-NH model; [26] for the large Inn valley in Alpine Austria with the COSMO model), though stable conditions often remain more difficult to simulate than convective boundary layers [27,28]. Our study is based on the Weather Research and Forecasting (WRF) model [29,30]. Starting from the version of previous research by [31] with a 1 km horizontal resolution, we developed a grid-nested mode version with the finest domain centered on the Cadarache site with a 111 m horizontal resolution. The model was run on a specific day corresponding to a well-documented case of the KASCADE-2017 campaign. The paper is organized as follows: The next section presents the site of the study and its geographical characteristics, the observations made during KASCADE-2017, the case study, and the numerical platform used for the simulation. The results of the simulation are then analyzed in line with the corresponding observations. The following section is a discussion about the performances of the simulation and the perspectives with which to improve it in the future. Our conclusions are presented in the final section.

2. Location, Observation, Case Study, and Numerical Platform

2.1. Location

The study area is located in the Alpine foothills of southeastern France. Figure 1 presents the five domains of the numerical simulation, which will be detailed in Section 2.4 and Section 2.5. This map allows these domains in western continental Europe to be situated. We will focus our study on the central part of the domain (enlarged in Figure 1b,c together with a representation of its topography). The main features are the mountainous areas of the Alps to the north and the east, the Mediterranean to the south, and the Rhône river valley to the west. We will focus especially on the Durance river (indicated as DV in Figure 1b), which flows to the SSW whilst separating plateaus of the pre-Alps. The main weather characteristics of this area are described in [24]. To summarize, the most frequent wind regimes are as follows: the “Mistral” wind, a cold and dry flow blowing from the north in the Rhône valley and then deviated according to the local topography. Southeastern flows, often moist and warm, bringing precipitation over the area. We then have spring-to-summer sea breezes, with an overall orientation from the south. These can be deviated along the valleys and deeply penetrate the Rhône valley as well as its tributary valleys [32]. The Mediterranean climate is characterized by frequent sunny and cloudless periods, with rainfall mainly occurring in the fall during short but sometimes heavy precipitation events.

At the local scale, which interests us (Figure 1c) the DV transitions from a ~6 km-wide, quite straight, and NNE–SSW-oriented valley (blue line in Figure 1b) to an E–W-oriented valley flowing towards the Rhône river, through an S-curve and a water gap (called the “Clue de Mirabeau”). Just upstream of this lies the smaller tributary, the Cadarache valley (CV), the primary focus of our simulation. It is SE–NW oriented (red line in Figure 1c), approximately 6 km long, 100 m deep, and 1–2 km wide, and its thalweg has a mean slope of 1.2°. The land use of the area is represented in Figure 1d, according to the Corine Landcover classification (https://sdi.eea.europa.eu/catalogue/srv9008075/api/records/5a5f43ca-1447-4ed0-b0a6-4bd2e17e4f4d, accessed on 10 August 2021). The CV floor and surroundings are mostly meadows with sparse trees (light green areas), non-residential buildings (purple), urbanized areas (red), and bodies of water (blue), whereas the hillsides and ridges are covered by denser vegetation from shrubs to sclerophyllous forests (medium and dark green).

2.2. Observations

The location of the various sites of observation is indicated in Figure 2. The first set of observations are continuous measurements routinely made at the low end of the Cadarache valley. The GBA station is a 110 m tower measuring pressure, relative humidity, and temperature at 2 m, as well as wind speed, wind direction, and temperature at 110 m, a level situated just above the adjacent ridgelines. In the central part of the CV, a Doppler sodar (Metek PCS.2000-24) is installed. This sodar measures the profile of the horizontal wind. The lowest height of observation is 60 m, and the highest is around 400 m, with a possible variation according to the propagation conditions of sound waves. At the same location, hereafter called “MET01” (see Figure 2), there is a weather station with observations at 2 m of atmospheric pressure, temperature, relative humidity, and incoming shortwave radiation, as well as wind speed and direction at 10 m.

From these continuous observations, it was possible to describe the overall cycle of the flows inside and just above the CV. Dupuy et al. [24] presented the wind roses at 10 m and 110 m, which revealed a frequent decoupling of the flows between the two levels: during stable conditions, SE winds follow the bottom of the CV; constrained between the valley sides in the first tens of meters in height, they are topped by NNE winds channeled by the larger Durance valley.

This superposition of flows on the vertical was observed in more detail during two field campaigns. The first one (KASCADE-2013) was conducted during the winter of 2012–2013. Thanks to regular profiling of the low atmosphere with radiosondes carried below either a tethered balloon or a free ascending balloon, it has been possible to build a consistent scheme of the complete cycle of these valley winds [21]. The second campaign (KASCADE-2017), on which we will focus in this paper, was conducted in the winter and spring of 2016–2017. The objective was to go beyond the first campaign and target horizontal variations in near-surface flows at the scale of the CV. For that purpose, a set of 12 weather stations was deployed during the entire campaign on both the bottom of the valley along its main axis, and on the flanks at various heights. To complement this, radiosondes were launched every three hours during the intense observation periods (IOPs) to document the evolution of the atmospheric profile. These local energy-budget measurement stations (LEMSs), developed at the University of Utah [33], measure wind speed and direction, incoming shortwave radiation, air temperature, and relative humidity at 2 m as well as atmospheric pressure at 1 m above the ground. In addition, surface temperature is measured with an infrared radiometer, and soil temperature and moisture content are measured at 5 and 25 cm depths. Furthermore, sonic anemometers were installed at four locations along the bottom of the valley and one in the DV, to document the turbulent structure of the flow close to the surface (2 m), through the record of the three wind components at a rate of 10 or 20 s⁻¹. Eventually, one site (AS1, close to GBA) was equipped with a four-component (short- and longwave, incoming and outgoing) radiation sensor. The location of these instruments is indicated in Figure 2. Eventually, a 2nd sodar (Remtech PA2) was installed close to the GBA tower, with the characteristics of 15 min, 25 m, 75 m, and 500 m regarding the time resolution, vertical resolution, and lowest and highest heights of observation, respectively. The KASCADE-2017 field campaign data are archived at: https://kascade.sedoo.fr, accessed on 10 August 2021.

2.3. Selection and Characteristics of the Case Study

Among the 14 IOPs which benefited from radiosonde profiles during KASCADE-2017, we selected IOP07, which lasted from 12:00 UTC February 20 to 12:00 UTC February 21. This choice was dictated by the availability of a complete set of observations during this period and suitable meteorological conditions for valley wind development during a stable stratification period (calm wind at the synoptic scale and cloud-free conditions). At a large scale, a high-pressure area was present over the Atlantic Ocean, and a low-pressure area over Central Europe, with a slight pressure difference between the two systems and therefore a very low pressure gradient over southern France (Figure 3). This W–E pressure gradient is favorable for the triggering of the Mistral wind in the Rhône valley, but with a moderate strength (average wind speeds of the order of 6–8 m s⁻¹ were observed, with gusts of up to 20 m s⁻¹), allowing the flow to enter tributary valleys of the Rhône river. During the daytime, the Mistral reaches locations as distant as our study area, where it arrived, in this case, abated and rotated to the west by the lowermost section of the Durance valley.

Figure 4 presents the time series of potential temperature, specific humidity, and wind, observed at a low level in the central part of the CV (MET01) and at the confluent zone of the DV and the CV (GBA). The diurnal cycle of temperature reflects a fair-weather situation, with a fast increase from dawn and a decrease initiated in mid-afternoon, accelerating just one hour before dusk. The diurnal temperature range (DTR) reaches around 18 K, which is only possible in cloudless, weak wind conditions. The wind and temperature observations allowed us to characterize the stability through estimates of a simplified bulk Richardson number, according to the formulation expressed in [23] (their Equation (1)). Such estimates return different results, according to the depth of the layer considered. The values computed at the GBA site, characteristic of the entire depth of the valley (110 m), reveal a very stable flow, with Ri values generally between 1 and 10 throughout the night, i.e., well above the critical threshold. On the contrary, the estimates at the MET01 site, computed in the layer below 10 m above ground, were below the critical threshold of 0.25, except at nightfall and during the period between 19 h and 23 h UTC. This means that the flow at the global scale of the valley is very stable throughout the night, while locally and near the surface, the stratification is closer to neutrality, and turbulent episodes can occur. We will discuss more about stability in Section 3.4. The wind speed always remains below 2.5 m s⁻¹ during IOP07. The day–night wind reversal is spectacular, with a westerly orientation during the daytime, behaving like a Mistral tail, as explained above, and an SE, CV-channeled katabatic wind during the stable stratification period [21], topped with a DV-channeled katabatic wind starting ~3 h after sunset. The specific humidity presents a complex evolution, with a rapid increase after sunrise. This could be related to the wind reversal occurring at this time, starting with vertical mixing, and/or to the start of evaporation with the availability of radiative energy. Other variations occur during periods of a stable wind direction and would require a more in-depth analysis in relation to the complexity of the wind field to be explained.

2.4. Configuration of the Numerical Simulation

2.4.1. Domain Settings

The WRF model [29] was run in order to simulate the IOP in a grid-nested mode. This mode allows simulations of both local and large scales of the flow while sparing computational resources and reducing numerical errors [34,35]. In this mode, the cell number of a given domain has the same order of magnitude whatever the domain. The cell size of the domain with the finest resolution is defined according to the scale of the local processes one wants to simulate, whereas that of the coarsest resolution is generally defined according to the resolution of large-scale forcing (operational analysis or reanalysis). The intermediate domains are staggered with a constant factor between the cell size of a domain and that of its parent and child domains. For WRF nested simulations, an odd factor has a better overlap with the parent domain, though this is not mandatory, as the two domains will have a common cell center [29], and lower factors require less interpolation.

WRF is routinely run in an operational forecast mode in southeast France by CEA with a horizontal resolution of 3 km. This resolution is too coarse to solve the thermally driven flows in the local valleys such as those observed during IOP07 [25]. Focusing on specific situations of KASCADE-2013, [31,36] have shown that WRF simulations with a kilometric horizontal resolution, though unable to resolve CV flows during stable conditions, are quite successful for the representation of the Durance down valley flow which develops in the lowest hundreds of meters above the ground some hours after sunset.

Since the CV, in which we are interested, has a width of 1–2 km, the horizontal resolution required to properly simulate the topographic forcing on the flows is O(100 m). For external boundary forcing, ECMWF reanalysis are available with a horizontal resolution of 27 km. We thus deployed five nested domains, named D1 to D5, with resolutions of 9000, 3000, 1000, 333, and 111 m, respectively. A constant factor of 3 is therefore kept between the resolution of two consecutive WRF domains, as well as between the coarsest WRF domain and the reanalysis. The 5 domains are displayed in Figure 1, the largest domain covering western continental Europe, and the finest one being an area of 18.5 × 17 km centered on the CV observations of KASCADE-2017.

For numerical constraints related to the two-way nesting mode, it is more practical to have identical vertical dimensions of the cells whatever the domain. The number and staggering of vertical levels are therefore defined to be consistent with the horizontal resolution of the finest domain, as well as with the physical processes we are planning to simulate. We defined 46 vertical levels, with stretching as one moves away from the surface. The lowest level center is at 6 m agl (above ground level), 9 levels are within the first 100 m agl, and 28 levels are within the first kilometer.

Table 1. Summary of the WRF domain and parameterization settings used for the simulations.

Date	2017-02-20 12:00:00 to 2017-02-21 12:00:00
WRF version	V 4.2
Global data forcing	ECMWF ERA5 1 h time step, 0.25° horizontal resolution, 38 vertical levels [37,38]
Nesting	Two-way
Vertical levels	46
Simulation time (h)	36
Spin-up (h)	12
Top of model (hPa)	50
Domain	D1	D2	D3	D4	D5
Horizontal resolution (m)	9000	3000	1000	333.333	111.111
Number of cells	106,100	100,100	100,121	175,178	169,154
Topography map resolution	5′	2′	30″	15″	3″
Time step (s)	45	15	5	0.5	0.125
Output interval (min)	180	180	10	10	10
Parameterizations
Microphysics	WRF single-moment 6-class scheme [39]
Planetary boundary layer	Quasi-normal scale elimination (QNSE) scheme [40]
Cumulus parametrization	Kain Fritsch [41]
Surface layer	QNSE surface layer unified [42]
Longwave radiation	RRTMG [43]
Shortwave radiation	RRTMG [43]
Land surface option	NOAH land surface model [44]

summarizes the model settings, where the time step varies according to the domain and is defined so as to satisfy CFL conditions. It ranges from 45 s for domain D1 to 0.125 s for domain D5. The model output is stored every 3 h for the D1 and D2 domains, and every 10 min for the D3 to D5 domains.

2.4.2. Representation of the Ground

We used maps with resolutions close to the resolutions of WRF and not too resource-intensive. From SRTM topographic maps [45,46], we extracted a 90 m resolution map for domain D5, whereas coarser resolutions were used for the other domains, up to 5 arcmin for the 9 km domain.

The representation of the land cover in the model has received particular attention, the methodology of which will be detailed in a forthcoming paper. Here, in order to go beyond the usual USGS land use dataset and its associated lookup tables, we took advantage of the Corine Land Cover (CLC) land use well suited for our European area, and we developed an aggregation method preventing misclassification when different land use categories are present in a single cell, at any given resolution. Although we developed new lookup tables (described in a forthcoming paper) assigning physical parameters (albedo, roughness length, etc.) appropriately suited to all CLC land use categories, in the present study, we used those derived for the USGS categories as devised by [47] so as not to burden the interpretations and to remain closer to existing studies using CLC.

Figure 5 illustrates how land use and topography are represented in the real situation, and in the 1 km and 111 m resolution domains. Starting from the greatest possible resolution on the left plots (the “ground truth” from the most detailed CLC map), we can see on the central plots that the 1 km resolution is unable to reproduce the real spatial heterogeneity both for land use and for orography. Furthermore, several weather stations deployed during the field campaign lie in the same cell; thus, the variability observed between these stations cannot be simulated. In contrast, at the finest resolution (right plots), the CV is clearly represented, and all stations are in different cells of the model.

2.5. Parameterizations

The choice of the schemes among the various possibilities offered by WRF mainly resulted from the sensitivity tests conducted by [31] in the simulation of a KASCADE-2013 case with a 1 km horizontal resolution. In substance, Table 1 mentions the settings that were adopted, and we provide a brief summary of our used parameterizations.

For land surface parameterization, we used the NOAH parameterization, which uses 4 soil layers and includes the freezing process, surface runoff, and a possible snow layer and calculates the soil moisture prognostically. We used the WRF single moment 6 for the microphysics; it prognostically calculates the mass of its six water species (rain, snow, water vapor, cloud water, cloud ice, and graupel). For both longwave and shortwave parameterization, we used the RRTMG model, which calculates the transmissivity of the atmosphere based on greenhouse gases and uses 14 wavelength bands. Lastly, for the PBL parameterization, we used the QNSE schemes, which include prognostic equations for the turbulent kinetic energy (TKE) and allow for TKE advection. They use a self-consistent, quasi-normal scale elimination algorithm and spectral space representation. We considered QNSE as a reference because previous studies have proven its good behavior at 1 km resolution [31], as well as at finer resolutions such as in the complex terrain of Sofia, Bulgaria, where QNSE has proven to be one of the most suitable parameterization schemes [48,49]. The boundary-layer parameterization merits further discussion, since it relates to the parameterization of subgrid-scale (SGS) turbulence. For kilometric (and coarser) resolutions, it is generally assumed that 3D turbulence occurs at scales which are lower than the horizontal dimension of the cell. The corresponding eddies are not solved by the model: the turbulence is entirely an SGS process and therefore is fully parameterized. At sub-kilometer resolutions, the largest eddies become of the same order as, or even greater than, the size of the cell, and a fraction of the TKE is thus explicitly solved by the simulation. The SGS turbulence scheme therefore has to be adapted, and when the major part of TKE is explicitly resolved, the model is run in the so-called large eddy simulation (LES) mode. Between fully parameterized and LES modes, there exists a range of scales belonging to the “gray zone”, where the dominant scales (i.e., the size of the most energetic eddies) are of the same order as the horizontal size of the cell, and the SGS parameterization is hard to define (Wyngaard [50] appropriately called this regime the “terra-incognita”).

We therefore attempted to estimate the size of the turbulent eddies in our case study. We know from the observations performed in the boundary layer during recent decades that the size of the energetic eddies increases with height, is larger in convective (i.e., daily) than in stable (i.e., nocturnal) conditions, and is commensurate with the boundary layer thickness [51]. Furthermore, flat and homogeneous areas are favorable for the development of coherent structures, which are of a size larger than that of structures over complex terrain such as our area.

We characterized the size of turbulence eddies through the estimate of the integral scale of the vertical wind computed from the signal recorded with a sonic anemometer installed at 10 m above the ground in the bottom of the CV (MET01 site in Figure 2). This value was then extrapolated up through the boundary layer thanks to the parameterization proposed by [52]. This parameterization requires the knowledge of the boundary layer depth (BLD), which was estimated from the radiosonde profiles. The method is explained, in detail, in Appendix A. In substance, during the daytime convective period, the BLD was of the order of 600–800 m, and the integral scales, computed from observations at 10 m, were in the range of 3–8 m and remained below 100 m when extrapolated through the boundary layer. Such results warranted the use of the fully parameterized SGS turbulence mode, even during the daytime part of the simulation. This choice was even more justified by the focus we have on stable stratification periods, where the characteristic scales are reduced with respect to daytime periods.

In order to bolster this choice, we also tested the SMS-3DTKE [53] scheme, released in WRF v4.2, which is a scale-adaptive TKE subgrid mixing scheme, extending the original TKE model [54] to the mesoscale and the gray zone resolutions. It functions, therefore, either as an LES or as a parameterized PBL scheme, depending on the domain resolution. In this test, the SMS-3DTKE scheme did not improve the wind at the lower levels of the CV or the global stability within the CV. Therefore, we will only discuss the results of the QNSE in the following sections. Similarly, the cumulus schemes (also known as convection schemes) have a gray zone from 8 km to 500 m. Weisman et al. [55] found that 4 km resolution grid sizes are fine enough for explicitly solving mesoscale weather. Jeworrek et al. [56] summarized that explicitly solved convection provides a better land–atmosphere interaction [57], while for coarser grids, parameterizations have better precipitation estimates. As we selected a case with neither rain nor clouds, we switched off the convection parameterization schemes at resolutions of 3 km and finer.

3. Results: Simulation vs. Observations

3.1. Improvement in the Refined Simulation

As mentioned in Section 2.3 and illustrated in Figure 4, IOP07 is characterized by a striking wind reversal, with a Cadarache down valley flow starting at dusk and persisting during the entire night. In spite of the light wind velocities (always below 2.5 m s⁻¹), the observed wind direction is quite steady. In order to evaluate the benefit brought by grid refinement, we performed three different runs of the model, in which the resolution of the innermost domain was either 1 km, 333 m, or 111 m. The results are illustrated in Figure 6 with the time series of the 10 m wind speed and direction at the MET01 location, both observed and simulated with various horizontal resolutions. All the results presented here are outputs of the finest domain of the run, preventing a two-way nested feedback to the parent domain. In the simulation with the coarse resolution (1 km), the model clearly fails to reproduce the wind direction during the night. This is not surprising and confirms the results obtained by [25,31,36]. Even with a 333 m resolution, the valley topography is too smooth, and thus the wind is not correctly simulated. In contrast, at the finest resolution (111 m), the characteristic dimensions of the valley, such as its depth, length, and width, are adequately represented, and the wind direction is correctly simulated, though the variability is higher than in the observations. During the daytime, the NW wind is well simulated in the afternoon of 20 February., but a discrepancy of 30 degrees–40 degrees is seen in the morning of 21 February. The performances of the simulations are characterized through the wind speed bias (with respect to the observations) and the direction accuracy (DACC45), the latter being the percentage of simulated directions which remain within ±45 degrees from the observation. These two metrics are computed over the 24 h period of the IOP. For the 111 m resolution, the DACC increases up to 60%, while the bias remains similar to the other resolutions around −0.32 m s⁻¹.

Another benefit from using a fine resolution is that the topography is close to reality, and a cell which contains a measurement station is in the model at an altitude (height above mean sea level; amsl) that does not differ from the true altitude of the station by more than 10 m, whereas at kilometric resolution, the topography is smoothed, and this difference can be as high as 50 m. Such altitude differences can have a considerable impact on the conditions of intense stratification, as observed during IOP07. For all these reasons, we will concentrate on the results of the finest (111 m) simulation in the following sections.

3.2. Overall Vertical Structure

Figure 7 presents the profiles of the potential temperature, specific humidity, and wind, both observed and simulated, at 18:00, 00:00 and 06:00 UTC during the IOP, roughly coinciding with sunset (17:13 UTC), the middle of the night, and sunrise (06:26 UTC). The radiosondes were launched from the MET01 location, and the simulated profiles were extracted from the finest domain simulation along the vertical in the cell that contains MET01. The profiles are represented below 1000 m amsl (i.e., 714 m agl at the site) only, since we are interested in the details of the low-level flows. Considering an ascending velocity of 5 m s⁻¹ for the balloon, and an average wind speed of 2 to 3 m s⁻¹ (as it was observed), the radiosonde observations presented in the figure are not strictly over the plumb line of MET01 but remain within a horizontal distance of 300–400 m, i.e., within the central part of the CV. It was found (not shown) that extracting the simulated values at the cells corresponding to the exact location of the balloon led to indistinguishable modifications on the profile.

The three profiles represent the early, middle, or late nighttime periods. Huge changes are observed from one profile to the next one, whatever the parameter. The 18:00 UTC potential temperature profile reveals that the surface layer is already strongly stratified, with an increase of about 6 K in the first hundred meters (observations), in spite of the short elapsed time since dusk. The surface continuously cools during the night, with a temperature 7 K cooler before dawn than after dusk, together with the boundary layer stratification, reaching 11 K between the surface and 300 m. At the top of the displayed profiles, we observe a cooling of the residual mixed layer, restricted to around 2 K during the night. Such cooling is not simulated by the model, which presents very close values at the top for the three profiles. With this restriction, the overall shape of the profiles is rather well reproduced. As in the observations, the model simulates a continuous cooling in the lower layers, together with an increasing stratification throughout the night. However, this cooling and the related intensity of the stratification are underpredicted, whatever the time. We will discuss low-level temperature simulation more in due course.

The moisture profiles are satisfactorily simulated at 00:00 and 06:00 UTC, with a surface layer which is, however, too moist (in excess by around 1 g kg⁻¹ with respect to the observations). In contrast, the 18:00 UTC profile is significantly too moist throughout the layer, exceeding the observations by 1–2 g kg⁻¹. This matter will be discussed in due course (Section 3.5).

The observed wind profiles reflect the development and evolution of the valley winds during stable conditions, as it was already described by [21]. At 18:00 UTC, the Cadarache down valley wind (direction in the range 115 degrees–155 degrees) has already started and persists throughout the night. It is restricted to a few tens of meters in depth, and its speed is below 2 m s⁻¹, as previously shown in Figure 4. NNE winds (Durance valley winds) reach the site around midnight, restricted below 200 m agl with weak velocities. The observations of the SOD1 (Figure 8), performed close to the GBA location, reveal that this flow strengthens and develops in height rapidly from this time, and during the 01–09 UTC period, the profile resembles that of 06:00 UTC, with a low-level jet (LLJ) of 4–5 m s⁻¹ in the 200–300 m agl layer. The development of valley winds is quite well simulated by the model, at least regarding the timing, wind direction, and vertical extensions (Figure 7 and Figure 8), even though the observations are at slightly higher altitudes than the model. The speed of the Durance LLJ is, however, underestimated on the 06:00 UTC profile (3 m s⁻¹ simulated instead of 5 m s⁻¹ observed), but examination of the sodar and model time–height cross-sections reveals that the agreement is rather good in the whole 00–06 UTC period. The major difference occurs at 18 UTC, with southerly winds simulated in the bulk of the 200–600 m agl layer (Figure 7), whereas westerly winds are observed there. In fact, the wind rotation from west to south does exist, as revealed by sodar observations. It is quite sudden and concerns the whole layer observed by the sodar (up to 500 m agl) but occurs around 19 UTC. The model simulates this rotation too early, which also explains the differences between observations and simulations in the temperature and moisture profiles of Figure 7. Another difference between observations and simulations is related to the deep layer of NNE wind (the Durance down valley wind, light blue color on the plot) occurring from around 20 UTC to the early morning of 21 February: it extends in the model down to heights lower than in the observations.

3.3. Diurnal Cycle

The diurnal cycle of the potential temperature and specific humidity at 2 m, as well as the wind speed and direction at 10 m, observed at the bottom of the CV (MET01 station—see Figure 2), is presented in Figure 9. While the wind evolution is rather well simulated, with an adequate wind speed and the reversal in the direction between the diurnal and nocturnal periods, we observe a large underestimation by the model of the daily temperature, whilst at night, it was overestimated. As a consequence, though the mean bias computed over the complete diurnal cycle is weak, the observed diurnal temperature range (DTR) is around 18 K, whereas only 12 K is simulated with a mean absolute error of 2.4 K. The main discrepancy occurs at the beginning of cooling: during the 3 h period starting 1 h before sunset, the observed temperature drops by around 12 K, whereas in the simulation, the corresponding cooling does not exceed 6–7 K. In contrast, during the rest of the night, the simulated temperature regularly cools down in a way fairly similar to the observations (−0.35 K h⁻¹). From sunrise on 21 February the simulated temperature increases less rapidly than observed and is therefore caught rapidly by the observed temperature, although heating rates remain similar from 10 UTC. DTR underestimation by WRF is a common problem and was already underlined in the same area by [25,31]. Lindvall and Svensson [58] even mentioned that DTR underestimation is a general problem in weather models. In our case, the origin of the cooling problem could lie in an incorrect representation of the soil properties, and/or a simulation of an air mass advection that differs from reality. We will speak more of stratification in a following section.

Concomitantly with the too weak cooling simulated in the early night, the model does not reproduce the evolution of the specific humidity. In the late afternoon, we observe a moistening of 2 h, followed by two hours of drying out, whereas the model continues to bring moisture to this location. The observed humidity always remains well below saturation, meaning dew deposition cannot be invoked to explain this decrease in the specific humidity. A possible explanation comes from an evening transition where turbulent vertical transport stops but evaporation continues, thus explaining the humidity peak. At first, humidity increases, since it is no longer transported away from the surface, whilst later, the humidity decreases when evaporation ceases or the nighttime regime starts, restarting some vertical transports [59,60]. A dip in the wind speed (Figure 9) and the reduced TKE (not shown) underpin a transport breakdown. For the remainder of the night, the observed and simulated moistures evolve similarly but maintain the difference created during the early night. From sunset on 21 February, as it is already observed in the temperature, the simulated moisture content increases less rapidly than in the observations, but after 10 UTC, the simulation and observations resemble each other more.

In the afternoon, wind directions match relatively well, although the observations are more spread out. The wind reversal of the evening transition occurs around the same time, after which the observations are constant around 135°, while the simulation shows more spreading, as previously mentioned. The timing of the morning transition does not match the observations precisely, as they maintain the valley orientation longer than the model. An earlier study showed that the wind direction change typically coincides with sunrise [21], resembling the current model behavior.

3.4. Simulation of the Stratification

For calm synoptic wind conditions, stratification is the leading mechanism for local wind development. In a vertical section crossing a valley, we could observe either horizontal isentropes, reflecting what is called a “mature cold pool”, or terrain-following isentropes, for “stratified terrain-following” situations, or an intermediate situation called “cold-air drainage” [61]. In the first situation, there are no slope flows; in the second, there are katabatic flows descending along the slopes whilst keeping the same potential temperature; and in the third situation, the downslope winds cross the isentropes and cool progressively. Mahrt [61] primarily characterized these different regimes through two parameters:

• The cold pool intensity (CPI), which is the potential temperature difference between the top of the flank and the bottom of the valley;
• The potential temperature difference on a horizontal plane between the top of the flank and the center of the valley (θ*).

For mature cold pools, θ* = 0, whereas CPI >> 0; for stratified terrain following, CPI = 0 and θ* << 0; and for cold-air drainage, θ* < 0 and CPI > 0.

In order to evaluate the performance of the model to simulate the stratification in the valley, we thus computed the θ* and CPI parameters from the observations and from the model. The CPI was estimated with the observations of two LEMS. One is situated at the bottom of the valley (B3, see Figure 2), while the other is situated on the northern flank, being either N3, which is situated at the top of the cross-valley instrumented axis, or N4, which is on the highest crest of this flank, but on an isolated hill and shifted to the north with respect to the instrumented axis. Meanwhile, θ* was estimated from the temperature difference between either N3 or N4 and the radiosonde values taken at the same altitude as the considered station. The observed θ* was therefore available every three hours, whereas CPI was observed every 5 min, since all LEMS have 5 min output intervals. CPI and θ* from the simulation were calculated at the same locations as for the observations and were available every 10 min.

The time series of the two parameters are presented in Figure 10, computed with either the N3 or N4 LEMS. The two estimates (with N3 or N4) are quite consistent. The observations reveal that the situation is close to that of a mature cold pool, with very high CPI values, reaching 7 K in the early night. The observed θ* values are close to zero or slightly positive, indicating a situation close to a “pure” mature cold pool (note that −θ* is represented instead of θ*, following the convention adopted by [61]). High values of the observed CPI prevent great cold-air drainage along the slopes of the valley flanks. On the contrary, the simulated CPIs increase around nightfall but remain much lower during the night, and θ* is weak or slightly negative during the night, but when combined, they reveal a possible drainage of the flow along the slopes of the valley flanks. We will discuss this question in due course.

3.5. Moisture Fields

In light of the model performance to simulate the winds and the stratification, we can now interpret the moisture fields. Figure 11 presents the horizontal fields of the simulated specific humidity close to the surface, on a horizontal domain of around 29 (lon.) × 25 km (lat.) surrounding the Cadarache valley, at 15:30 and 18:00 UTC, i.e., at the end of the diurnal convective period and at the beginning of the stable stratification period. The fields of the simulated wind direction at 10 m are also drawn at the same times. In the late afternoon, the E to SE winds establish along the valleys/ridges in the central part of the represented domain (see the purple and violet bands on the wind direction 18:00 plot). As a consequence, the area with higher specific humidity observed in the SE corner at 15:30 is advected throughout the domain to the NW and reaches the Cadarache valley. We do not discern any variation related to the local topography at the 15:30 fields; in other words, the values are about the same at the bottom of the local valleys or at the neighboring crests.

This is also illustrated in the time series simulated throughout the IOP, which is very comparable at the bottom of the valley or at the top of its flank (Figure 11g,h). In contrast, the observed values at the same locations reveal a considerable difference: while the plot at N3 resembles the simulation output, the B3 values at the valley bottom are significantly lower during the night. One hour before sunset, the observed moisture decreases at the bottom of the CV, which is not simulated by the model. Meanwhile, on the crest, we also observe a decrease starting at the same time as that at the bottom of the valley, but this is rapidly compensated for by a sudden increase occurring one hour after sunset. At this location, the model matches the observed time series quite well but anticipates the sudden increase by around 1.5 h.

The differences between the observed and simulated time series at the bottom of the valley can be related to the differences between the observed and simulated stratification in the valley. As mentioned above (Section 4.4), the observations reveal quite a mature cold pool situation, which implies that there are almost no downslope flows along the flanks of the CV. The air mass observed at MET01 during the night, therefore, was constrained close to the bottom of the valley where it circulated following the valley’s orientation. In contrast, the model simulates a mix of such a flow with downslope currents along the valley flanks, which tends to homogenize the air physical properties at the various levels.

The timing of the humid air arrival seems to affect the cold pool development. Comparing Figure 10 with Figure 11h shows that the CPI increases until the arrival of the moist air mass. In the simulation, the CPI growth stops earlier, coinciding with the early humidity arrival (at 17 UTC), and the simulated CPI remains constant thereafter. The observations in Figure 10, however, show that the observed CPI continues to increase, up to the observed arrival of the humid air mass. At that moment (a little before 19 UTC), the CPI shows a dip but restores relatively quickly and does not increase further. In both the simulation and observation, the humidity arrival matches the moment with the largest CPI values and occurs earlier in the simulation.

3.6. Vertical Cross-Sections

Three vertical cross-sections of the wind speed and direction at 00 UTC are presented in Figure 12, along the tracks indicated in Figure 1 and Figure 11. The top one, 8 km long along the Durance valley, shows a down valley flow, 300–400 m thick, with a well-marked LLJ of around 7 m s⁻¹ at 150–200 m above the ground. The shear layer, between the down valley wind and the regional-scale westerly flow lying above, is quite thin and slightly rises in the up-valley direction. The apparently south-originating flow layer at 700 m amsl in Figure 12b, and d is an artifact originating from the interpolation, which cannot handle cyclic color scales and will linearly rotate from NNE to NNW (0° to 359°) flows. For this reason, the simulated wind at 5000 to 7000 m along the CV cross-section (Figure 12d) shows the slow transition from NNE to NNW, with a sharp artifact change at the 0° to 359° transition. The terrain elevation, close to the lower end of the DV section, corresponds to a protrusion of the S-shaped right bank of the valley, with a cliff-like sidewall about 100 m high (Figure 1). This wall guides the DV katabatic flow and forces it to follow the river bend. Similarly, on the downstream side of this obstacle, the patch of low-altitude SE flow reflects the direction of the S-shaped valley, where it turns to the northwest (close to the Clue de Mirabeau), possibly reinforced by another down valley flow from the small valley parallel to the CV. At higher levels (i.e., above 400 m asl), the flow is hardly disturbed by this obstacle, as shown by the continuation of the LLJ which is only ascending marginally. It is interesting to note that the CV flow is not visible here (the CV axis crosses the section at an abscissa of about 2 km, indicated by the red arrow), meaning that this flow does not penetrate the Durance valley at that point. Upstream in the DV (at abscissas between 4 and 6 km) and very close to the surface, we can observe weak winds from the west, probably resulting from down valley flows in the many small tributary valleys discharging on the right bank of the DV.

The section along the CV confirms that the down valley flow is weak, confined to a shallow layer (a few tens of meters), and does not reach the DV at this time. The Durance down valley flow is present all along the CV, just above it, but the LLJ is restricted to a zone of around 3 km wide, centered above the Durance river. In the section perpendicular to the CV axis (Figure 12e,f), we see that the Cadarache down valley flow is mainly present in the central part of the valley, where its depth reaches around 100 m. The figure also reveals the asymmetrical behavior of the flow on the two sidewalls, with a thicker (~100 m) layer of very weak winds on the southern flank than on the northern flank. With such an orientation, these winds in the model descend the grade and contribute to supplying the thalweg of the CV with air masses from above, as we suggested earlier.

4. Discussion

4.1. Wind Modeling—Vertical Profiles

With respect to the previous WRF simulations conducted in this area with a kilometric resolution, be it operational simulations as analyzed by [25], or case studies such as in [31,36], we demonstrated that a simulation with a hectometric resolution is clearly able to reproduce, with some accuracy, the observed winds and their diurnal cycle in the small-scale valley, which was unreachable with coarser resolutions. The main driving mechanism of the local flow, i.e., topographic forcing, is adequately simulated by the model. As a consequence, the diurnal cycle of the wind profile, as observed by a sodar, and the shear layer, during nocturnal stable periods, between the Cadarache down valley SE flow, close to the surface, and the Durance down valley NNE flow above, are well simulated. We can now consider that a reliable tool is available for the simulation of the vertical structure of the flow as it has been described by [21] thanks to the observations of the KASCADE-2013 field campaign. This WRF platform extends the predictive capability of already available tools, be it the simple diagnosis of the Cadarache down valley wind of [23], the nowcasting based on downscaling from routine observations of [24], or the forecasting based on operational coarse simulations and ANNs by [25].

Vertical cross-sections show a well-developed Durance down valley flow, stretching over the CV, which matches the observations of [21], even though the simulated Cadarache flow is shallower than the observations. Regarding the Cadarache down valley flow, the simulations show a smaller layer on the south flank than on the north flank, which could not have been observed with observations on a single vertical.

The Durance LLJ has limited influence on the winds close to the surface, in contrast to what was reported by [62] in relatively uncomplicated terrain (a valley leading to a plateau in the Pyrenean foothills). Nevertheless, further up the valley, the jet in [62] appears to break away from the surface near the confluent zones of wider valleys, which resembles the situation around the confluent zone of the CV and DV. Such jet streams are sometimes associated with internal waves, as reported by [63] through Doppler lidar observations, but such events were beyond the capabilities of our observations.

4.2. Stratification and Slope Flow

The sensor network deployed at the surface during the KASCADE-2017 campaign, together with the profiles obtained from radiosondes, allowed us to characterize the stratification inside the valley. We must admit that the model does not reproduce the temperature field structure accurately enough. In line with the underestimation of the DTR, already highlighted by [31], the cold pool predicted during the night is far from being as mature as that revealed by the observations, agreeing with the findings of [28] on the difficulties of creating a cold pool in fine-scale simulations and demonstrated with the values of the CPI and θ* parameters. An effort had been undertaken to improve the representation of the land use in the model, but the consequences on the DTR remained modest. As a consequence, the simulated temperature differences between the bottom of the valley and the top of the sidewalls (CPI) were underestimated, with the model generating downslope flows which are probably not realistic given the observed CPI values. We have to note, however, that the observations on the sidewalls were insufficient for a full diagnosis of downslope flows, since (1) only horizontal wind components were measured, whereas a complete diagnosis would require 3D anemometers at the locations of the LEMS stations; and (2) the wind velocities observed during the nighttime periods were often very weak, i.e., below the threshold for which the cup and vane anemometer measurements can be considered as reliable. For future improvements, we are planning to incorporate into the model a better representation of the soil characteristics. In particular, the soil water content is known to play a crucial role in the surface energy partitioning (evaporative fraction), and previous studies have shown that ECMWF analyses overestimate soil moisture, especially for dry areas [63]. The consequence is a too low surface sensible heat flux and surface temperature during the daytime.

4.3. Moisture Field and Transport

To summarize, the horizontal advection is quite well simulated, whereas the vertical transport (along the slopes) is not. As a consequence, the specific humidity is vertically more homogenized in the simulations than in reality. The air mass observed at a given station at the bottom of the valley comes from an upstream area close to the ground, whereas the model mixes this air with downslope currents. The observed humidity is therefore much more dependent on the local terrain over which the down valley flow has passed.

4.4. Perspective

The perspectives of this work are therefore to improve the model performance in the representation of the DTR and the stratification. In particular, we will explore the effect of the moisture in the soil, which could be overestimated in the simulation, leading to a reduced DTR. Once we are confident enough in the performance of the model, the goal will be to generate a high enough number of fine-scale simulations to train ANNs on as many locations as necessary to downscale the 1 km operational weather simulations and generate a reliable representation of the flow in a complex area whilst avoiding a computationally expensive operational, fine-scale simulation.

Kalverla et al. [31] found that the QNSE is a good PBL parameterization for this location at 1 km resolution grids. Sensitivity tests with different PBL schemes at sub-kilometer resolution under stable conditions benefited only from limited research because most studies focused on the daytime simulations, where the turbulent scales are larger and LESs are required. PBL parameterization at these resolutions in stable conditions remains a very open question.

5. Conclusions

The study presented in this paper is a contribution to an overall approach aiming at the development and qualification of various tools intended to diagnose and/or forecast the flows at the very local scale in complex terrain areas. It focuses on stable stratification conditions because they are often more penalizing regarding the consequences of chronic or accidental release of pollutants or contaminants, due to the fact that vertical dilution by turbulence is reduced. Furthermore, combined channeling and stratification effects in complex terrain challenge the monitoring of plume behavior and render conventional transport and dispersion models based on crude assumptions regarding, e.g., horizontal homogeneity and wind field structure ineffective.

We thus developed, over recent years, various tools, ranging from very simple diagnostics of valley flows based on routine observations [23], or statistical methods based on the ANN technique applied either to routine observations [24] or operational forecasts [25]. The ANN technique is certainly very efficient in nowcasting/forecasting the wind at a given spot; however, it remains difficult to extend to a large number of locations because the learning process would have to be based on a long series of observations at each targeted location (“long” meaning here that all possible weather conditions have been encompassed during observation periods).

For these reasons, we engaged in the qualification of a high-resolution weather model able to solve the flows at the very local scale. This model is based on the WRF platform used worldwide, which is already run in an operational mode at a kilometric scale resolution in southeastern France. The version developed during this work has five nested domains, the horizontal resolutions of which range from 9 km for the coarsest to 111 m for the finest. The latter is centered on the area we are interested in, the central part of which has been widely documented during the two field campaigns of KASCADE-2013 and KASCADE-2017. These observations revealed that, during stable stratification periods with calm to moderate winds, the low-level flow has a well-marked diurnal cycle, with, in particular, a shallow down valley wind layer starting at sunset in the small Cadarache valley.

The model was run on the IOP07 of the KASCADE-2017 campaign. The 24 h simulation encompassed the afternoon, the entire night, and the following morning. When confronted with the observations, the model revealed its capacity to simulate the diurnal cycle of the flow. In particular, the nocturnal down valley flow in the 1 km-wide CV is well reproduced, regarding its strength, orientation, thickness, and diurnal cycle. This constitutes considerable progress with regard to the performances of previous WRF simulations performed in this area, be it operational weather forecasts or dedicated studies conducted on specific days with kilometric resolution. The good comparison with sodar profiles shows that the flow is also well simulated above the CV, thus indicating that the simulation is also relevant at scales greater than those of the CV.

As it stands, the model, however, suffers from some weaknesses. The diurnal temperature range is underestimated, and as a consequence, the stable stratification does not have the appropriate intensity. Rather than a mature cold pool in the CV, as suggested by the observations, the model simulates a mixed cold pool, with significant drainage along the sidewalls, whereas the observed winds do not show slope flows, at least in the moments with wind speeds high enough to activate the cup anemometers. As a consequence, the model overestimates the transport of scalars, such as the specific humidity, from the sidewalls toward the bottom of the valley, which is not observed.