Boundary‐layer plumes over mountainous terrain in idealized large‐eddy simulations

Coherent plume structures in the convective boundary layer over non‐flat terrain are investigated using large‐eddy simulation. A conditional sampling method based on the concentration of a decaying passive tracer is implemented in order to identify the boundary‐layer plumes objectively. Conditional sampling allows quantification of the contribution of plume structures to the vertical transport of heat and moisture. A first set of simulations analyzes the flow over an idealized valley, where the terrain elevation only varies along one horizontal coordinate axis. In this case, vertical transport by coherent structures is the dominant contribution to the turbulent components of both heat and moisture flux. It is comparable in magnitude to the advective transport by the mean slope‐wind circulation, although it is more important for heat than for moisture transport. A second set of simulations considers flow over terrain with a complex texture, drawn from an actual digital elevation model. In this case, conditional sampling is carried out by using a simple domain‐decomposition approach. We demonstrate that thermal updrafts are generally more frequent on hill tops than over the surroundings, but they are less persistent on the windward sides when large‐scale winds are present in the free atmosphere. Large‐scale upper level winds tend to reduce the vertical moisture transport by the slope winds.

turbulent fluxes due to these flow structures are not explicitly resolved and need to be parametrized.Thereby, the interest in coherent convective structures was often driven from a cloud physics point of view (Simpson and Wiggert, 1969;Nicholls, 1989).Generally, in parametrizations of the CBL for weather and climate models, the traditional eddy-diffusivity formulation is only able to model local turbulent fluxes (Siebesma and Cuijpers, 1995).The properties of convective updrafts in the daytime boundary layer are often poorly correlated with local gradients; hence, they need to be parametrized using non-local turbulence closures (Randall et al., 1992).To this end, the eddy-diffusivity approach can be combined with a mass-flux term to the so-called eddy-diffusivity mass-flux scheme providing a better physical basis for parametrizations (Siebesma et al., 2007;Angevine et al., 2010).
Unlike NWP models, large-eddy simulation (LES) can be used to explicitly resolve small-scale coherent motions down to the meter scale.They have been widely used to study both clear and cloudy boundary layers, often with the goal to improve parametrizations (Schumann and Moeng, 1991;Heus and Jonker, 2008) and to investigate the spatial characteristics of dry convection (Babi ć and De Wekker, 2019).Statistics of LES output can reveal how the total turbulent flux is partitioned between local (diffusive) and non-local (related to coherent structures) components, provided that a method for the accurate identification of coherent structures is available (Raupach, 1981).Several updraft sampling methods relying on vertical velocity and thermodynamic variables have been introduced and applied in dry CBL cases (Williams and Hacker, 1992;Berg and Stull, 2004).Owing to the limitations of such approaches in the transition between sub-cloud and cloud layer, sampling methods based on passive tracers were introduced, with the aim to evaluate the entrainment and detrainment rates into the plumes from the surface to the top of the cumulus clouds (Couvreux et al., 2010).Other publications point out the importance of taking into account downdraft structures as well; for example, in stratocumulus-topped marine cases (Brient et al., 2019).It has also been shown that coherent convective structures in simulation data can be analyzed with the help of joint probability density functions (JPDFs; Chinita et al., 2018).The JPDFs of the vertical velocity with temperature and humidity perturbations have often been used to study and to parametrize dry and especially moist convection (Wyngaard and Moeng, 1992;Larson and Golaz, 2005).Non-local, coherent turbulent fluctuations manifest in the non-Gaussian part of the JPDF.
Most existing studies about coherent plume structures are limited to the flat boundary layer.The term "mountain boundary layer" refers to the boundary layer over complex terrain responding to surface and terrain forcings with time-scales of a few hours responsible for the exchange of heat, mass, and momentum between the mountainous terrain and the free atmosphere.In contrast to flat terrain, this exchange results from multiple processes on different spatial and temporal scales, such as orographic gravity waves, turbulent motions, moist convection, and thermally induced circulations.This implies substantial modeling challenges due to lack of parametrizations and the need of high grid resolutions (Lehner and Rotach, 2018).Based on high-resolution model simulations, the reconstruction of transport processes in a valley is possible after a spatial and temporal filtering of the flow allowing a decomposition into a large-scale wind, a local mean circulation, and a small-scale turbulent part.This approach has led to numerous insights related to the advective transport by local thermal wind systems, such as slope and valley winds, in particular during daytime.In the case of the slope-wind circulation, a buoyancy-driven upslope flow arises over the heated valley slope (Vergeiner and Dreiseitl, 1987).A subsiding branch running from the ridge top back to the valley center closes the circulation cell.Budget analyses show that the slope-wind circulation has a cooling effect over the slopes and it leads to local heating in the valley center (Schmidli, 2013).Together with the turbulent heat flux, this results in a rather homogeneous heating of the valley atmosphere.Moisture, approximated as a passive scalar quantity, is efficiently exported out of the valley by the slope winds and accumulated over the ridges (Kuwagata and Kimura, 1997;Weinkaemmerer et al., 2022).In consequence, it is found that a large-scale, horizontal wind effectively reduces the export of moisture out of the valley by mixing dry air from above the valley center into the moist slope-wind layer.In contrast, the temperature distribution remains nearly unaffected.Unlike a horizontally homogeneous CBL, where buoyant plumes are randomly distributed and feature a broad spectrum of scales (Griewank et al., 2022), the CBL over orography can be expected to show particularly intense plumes over the hilltops, which constitute natural regions of horizontal convergence in daytime.Also, their spatial scales may be controlled by the properties of the underlying orography (Göbel et al., 2023).Several studies focus on the vertical transport induced by stationary updrafts over mountain crests, primarily because they are the most likely to evolve into precipitating moist convection (Demko and Geerts, 2010a;Demko and Geerts, 2010b;Nelson et al., 2022;Göbel et al., 2023).
Owing to an upper level wind, they can be shifted horizontally (Panosetti et al., 2016).Aside from the hilltop updrafts, it is not yet clear how the boundary-layer plumes behave in the slope-wind layer and how the turbulent fluxes are partitioned between coherent-motion transport and unorganized turbulence at different locations along the valley cross-section.
In this article, we want to address these unclarities by studying boundary-layer plumes over idealized and semi-idealized complex orography.For this purpose, we run LES at about 40 m resolution and implement a conditional sampling method that relies on a transported passive tracer.Following Weinkaemmerer et al. (2022), basic experiments are performed over a quasi-two-dimensional, periodic valley where the surface fluxes are prescribed and constant.This idealized-valley case simplifies the updraft sampling and the decomposition of the flow.The goal is to identify and characterize boundary-layer plumes over smooth complex terrain, to study the horizontal distribution of the turbulent flux contributions, and to gain insight into the structure of the slope-wind layer.Furthermore, the turbulence statistics in relation to the local mean circulation are evaluated by calculating the JPDFs on different height levels.This is followed by two simulations over three-dimensional orography with a radiation-driven diurnal cycle.Thereby, a simple approach towards structure identification and flow decomposition over more realistic terrain is tested.With that, the effect of an upper level wind on the coherent structures and on the bulk quantities is analyzed.
The article is structured as follows: Section 2 presents the experimental set-up and the numerical model.The analysis methods with the conditional sampling and the flow decomposition are introduced in Section 3. Sections 4 and 5 discuss the results of the simulations.Conclusions are given in Section 6.

Experimental set-up
The idealized-valley case is set up as in Weinkaemmerer et al. (2022).The orography is a periodic, infinitely long sinusoidal valley described by where z s is the surface height, h = 1,500 mm is the ridge height from the valley floor to the crest, and W = 20.48km is the width of the valley from ridge to ridge, which also equals the domain width L x .The domain length L y is 40.96km.In the reference experiment (REF), the background atmosphere is at rest.The initial atmosphere is in hydrostatic balance with p = 1,000 hPa and  = 297 K at z = 0, a linear potential-temperature increase with height of Γ = 0.003 K⋅m −1 , and a constant relative humidity of 40%.The surface sensible heat and moisture fluxes are respectively prescribed as (w ′  ′ ) s = 0.12 K⋅m⋅s −1 and (w ′ q ′ ) s = 0.05 g⋅kg −1 ⋅m⋅s −1 (corresponding to a Bowen ratio ≈1), and the roughness length is set to z 0 = 0.16 m.The total duration of the simulation is 6 hr.An additional simulation is carried out over a flat surface with the same set-up as is used over the valley.
For the three-dimensional cases, the terrain is loosely derived from the area around the Chabre mountain near Laragne-Montéglin in southern France, standing as an exemplar for a hilly terrain.The smoothed orography is damped towards a constant height of 600 m above sea level at the lateral boundaries of the domain, where periodic boundary conditions are imposed.The dimensions of the domain are L x = 34.9km and L y = 32.3tkm.The total duration of these simulations is 12 hr.The initial sounding corresponds to July 11, 2003, 0600 UTC, a summer day featuring a clear CBL.For the three-dimensional reference case (REF3D), the winds are set to zero.For the WIND3D case, the full initial sounding is adopted, which includes zonal and meridional wind.The time and the location determine the surface fluxes given by a surface-layer model in combination with a radiation scheme.A brief overview over the different simulations is given in Table 1.

Numerical model
The numerical simulations are performed with Cloud Model 1 (Bryan and Fritsch, 2002) in LES mode.It uses terrain-following  z coordinates.The model is integrated using third-order Runge-Kutta time differencing (Wicker and Skamarock, 2002) and a fifth-order weighted essentially non-oscillatory scheme both for the advection of scalars and momentum (Jiang and Shu, 1996).The Coriolis force is neglected.The horizontal resolution is 40 m (42 m for the three-dimensional cases).Grid stretching is applied in the vertical, reaching from 8 m vertical resolution at the surface to 40 m at a height of 2,800 m and above.The model top at 5,000 m is a rigid lid.This corresponds to 175 vertical levels.For the WIND3D case, the model top is raised to 8,000 m (250 vertical levels) because of vertical gravity-wave propagation.In order to suppress spurious wave reflections, a sponge layer with Rayleigh damping is employed ( = 300 s).Its depth is 1,000 m (2,000 m for WIND3D).In order to facilitate the onset of convective turbulence, the initial atmosphere is altered by a random temperature perturbation with a maximum amplitude of 0.5 K.The model uses an adaptive time step with a typical value of 2 s.For subgrid-scale turbulence closure, a 1.5-order turbulent kinetic energy scheme is chosen (Deardorff, 1980).Microphysical processes are parametrized using the Morrison double-moment scheme (Morrison et al., 2009).Atmospheric radiation is parametrized by the Goddard scheme (Chou and Suarez, 1999).For the surface-layer model (Jiménez et al., 2012), the grassland option is selected.

Conditional sampling
The conditional sampling method has been implemented according to Couvreux et al. (2010).It uses a passive tracer emitted at the surface with a constant surface flux.The value of the surface flux is arbitrary.The tracer concentration decays exponentially with a time-scale of  0 = 15 min to ensure that the tracer does not accumulate outside the updrafts above the surface.A grid point is characterized as belonging to a coherent structure if the local vertical velocity w is positive and the tracer concentration anomaly c ′ exceeds the standard deviation of the tracer concentration  c and a minimum threshold  min : In Couvreux et al. (2010), c ′ is defined as the deviation from the domain average of the tracer concentration.This is only meaningful in a horizontally homogeneous set-up; that is, over flat terrain.In this study, for the idealized periodic valley, the spatial average is computed along the y-axis of the domain; that is, the symmetry axis of the terrain.Sensitivity tests have shown that a number of grid points of ≈1,000 is sufficient for the spatial average to converge.With that, the conditional sampling is applied every 1 min of the simulation.The minimum threshold is defined as selected in a non-turbulent environment.The effect of an additional condition, ensuring that only cloudy grid points are selected within the cloud layer, is negligible in this study as only sparse cumulus clouds are formed in the REF case after 4 hr simulation time.
Over three-dimensional orography, the computation of a spatial average is not straightforward without a symmetry axis of the terrain.Thus, the domain needs to be divided into multiple subdomains (see also Rai et al., 2017a;Rai et al., 2017b;Babić and De Wekker, 2019).The spatial (horizontal) averaging is performed on the  coordinate surfaces, over non-overlapping subdomains of size 1.3 × 1.3 km 2 .This approach is a compromise between the competing needs of obtaining robust statistics of coherent structures (requiring large subdomains) and considering areas of relatively homogeneous orography (requiring small subdomains).It ensures that the averaging still extends over 1,024 grid points.In terms of its size, each square subdomain can be considered as representative of a grid element of an operational-scale NWP model.

Flow decomposition
By means of spatial and temporal averaging, the flow is decomposed as shown in Weinkaemmerer et al. (2022).
With the aim to compute the turbulent fluxes, the turbulent fluctuation a ′ of a flow variable a is separated from the ensemble average a via Reynolds decomposition: (2) Subgrid-scale fluctuations not resolved in the LES are treated as being included in the turbulent part.Following the Reynolds averaging rules, the ensemble average of a product is given by The second-order turbulent fluxes consist of a resolved and subgrid part: For the analysis, the subgrid parts are obtained from the LES subgrid turbulence model.
In order to study the thermally induced flow, the quantity's horizontal average over the whole domain ⟨a⟩ can be removed from the ensemble average to leave a local average a c describing the mean circulation: (5) Over the sinusoidal valley, the required ensemble average is approximated by an average in time and in the along-valley direction: where x, y, and z are respectively the eastward, northward, and vertical directions, and t is time.Here, z denotes the absolute height above sea level.The time averaging is based on 1 min intervals.It is done online during the model run in order to reduce the data output.Each averaging period T extends over 40 min as a compromise between accuracy and stationarity (Schmidli, 2013).The domain mean is then achieved by averaging a additionally in the x direction on constant height levels: This requires interpolation in the z direction.Over three-dimensional terrain, the flow components are calculated analogously.For the ensemble mean, the averaging is performed on the  levels over each subdomain: where A sub is the subdomain area.Then, the domain average over the total domain area A tot gives a vertical profile again:

Joint probability density functions
The vertical turbulent flux of a variable  is determined by the covariance of its turbulent part  ′ with the perturbation of the vertical wind w ′ .Here, we refer to the turbulent flux averaged over the x-y plane at a given height z above sea level.In the case of discrete model data, this can be expressed by the sum over all N data points at that height: The use of terrain-following coordinates implies that model levels are not horizontal, so vertical interpolation is required to obtain w ′ and  ′ from the LES data.A continuous probability density function can be found to describe the two-dimensional histogram of w ′ and  ′ so that where f (z, w ′ ,  ′ ) is the so-called JPDF, which is normalized to unity.In this study, the JPDFs are fitted to the model data via a Gaussian kernel density estimation.With the help of the updraft sampling, a JPDF for the contribution of the plumes to the total turbulent flux can be specified as well.Figure 1 shows the turbulence statistics for the turbulent heat and moisture flux over flat terrain for demonstration.Per definition, the JPDF for the non-local mixing by the coherent structures (red) corresponds to positive vertical-wind components.The small but still positive densities for w < 0 are an artifact and a drawback of the Gaussian kernel density estimation.At the selected time, the boundary-layer height is about 1,500 m.Especially in the lower half of the boundary layer, the probability density is highly skewed towards the first quadrant, meaning that positive potential-temperature perturbations are transported upwards by the plumes.At 800 m (Figure 1b), the coherent structures are responsible for almost 90% of the turbulent heat flux.At higher altitudes, the overshooting thermals are potentially cooler than their surrounding environment.Thus, their contributions to the heat flux are located in the fourth quadrant (Figure 1d).The remaining second and third quadrants comprise contributions from the turbulence outside the plumes; that is, downward turbulent motions resulting in negative (positive) flux contributions when they are correlated with positive (negative) potential-temperature perturbations.Thus, the second quadrant is especially important in the entrainment zone (Figure 1d), and positive contributions from downdrafts are generally located in the third quadrant.For the turbulent moisture flux (Figure 1e-h), the JPDF for the coherent transport is generally located in the first quadrant, as moisture is transported away from the moist surface layer to the drier air above.Thus, the moisture-enriched coherent structures are more dominant at higher levels.

Boundary-layer plume characteristics
In this section, the results of the conditional sampling will be characterized for the idealized-valley case.As an example, Figure 2a shows a cross-section of the valley at an arbitrary point in the y direction.This is , , after a spin-up time of 4 hr when the simulation is expected to represent typical midday conditions.Following Schmidli (2013), we distinguish between mixed-layer height and boundary-layer height.The height of the nearly neutral turbulent mixed layer is obtained as the height where the vertical gradient of the along-valley averaged potential temperature first reaches a critical value of  c = 0.001 km −1 .The mixed layer is topped by a stable entrainment zone.The top of the entrainment zone equals the height of the atmospheric boundary layer, which is diagnosed via a simple maximum-gradient method as in Schmidli (2013).
It can be seen that plume structures have evolved all along the valley cross-section, whereby their vertical extent scales with the boundary-layer height.The tallest structures can be found in the valley center and over the ridges.As the boundary-layer height is a spatially averaged quantity, overshooting plumes also occur.To some degree, elevated structures that are decoupled from the surface also appear above the boundary-layer height.A comparison with Figure 2b shows that these structures are located in the recirculation zone above the ridges.Figure 2b also reveals that the slope-wind circulation has formed two vertically stacked circulation cells, as is typical for deep valleys (Wagner et al., 2015).Roughly, the mixed-layer height coincides with the height of the upslope branch of the cross-valley circulation.The highest horizontal wind velocities prevail at higher altitudes in the slope-wind layer close to the ridge.In top view, Figure 2c confirms that the deepest structures can in fact be found over the valley bottom, where the boundary layer is deepest and the slope winds are weak.Qualitatively, the boundary layer in the valley resembles most closely the horizontally homogeneous boundary layer over flat terrain.
It is a goal of this study to gain insight into the characteristics of the plumes in the slope-wind layer and how they interact with the slope wind.From one time step to the next, the plumes seem to move upslope while they are evolving and decaying.Instead of identifying individual convective objects, as in Griewank et al. (2022), and then tracking their motion, we estimate the bulk velocity of coherent structures with an image correlation analysis (Figure 2d).Thereby, the sampled structures are shifted in slope direction until the spatial correlation with the structure pattern of the following sampling step (1 min ahead) reaches a maximum.This is only done in the lowest vertical model layer and in intervals of 1 km along the x-axis.From the extent of this horizontal shift, a bulk velocity of the coherent structures can be calculated.The result shows that the horizontal velocity of the coherent structures over the valley slope is comparable to the surface wind speed.At higher altitudes, this technique does not produce reliable results because of the rapid changes in the shape of the structures due to turbulent motions.However, as the plumes between x = 4 km and x = 7 km are mostly limited to the slope-wind layer, it can be assumed that their bulk transfer velocity does not differ substantially from u c .

Heat and moisture flux decomposition
Based on the conditional sampling, the contribution of the coherent structures to the vertical turbulent flux of heat and moisture can be quantified.To this end, the turbulent flux is decomposed into a coherent-structure part (top-hat representation), a second term expressing the in-structure variability of the updrafts, and a third term representing the turbulence in the environment outside the plumes (Siebesma and Cuijpers, 1995;Couvreux et al., 2010): where  is the local coherent-structure coverage, superscript "i" stands for in-structure, and superscript "e" stands for environment.This means that the calculation of the ensemble mean according to Equation ( 6) is carried out either over the grid points inside or outside the structures.The resulting cross-section plots are shown in Figure 3.It is striking that the in-structure variability of the turbulent fluxes is almost negligible both for the heat and the moisture flux.This confirms the coherency of the sampled updrafts.Generally, the total turbulent flux (Figure 3a,e) is dominated by the plumes.For the heat flux, Figure 3b shows that the positively buoyant coherent structures generally cause a positive turbulent heat flux in the lower half of the boundary layer all along the valley cross-section.In the entrainment zone, the heat flux turns negative due to overshooting plumes that are cooler than their environment in terms of the potential temperature.The environmental turbulent heat flux counteracts the effect of the plumes to some extent, except in the lower boundary layer, where it has a positive contribution as well (Figure 3c).At the boundary-layer top, for instance, this traces back to descending positive potential-temperature perturbations and ascending negative ones.A different picture is shown by the turbulent moisture flux, where the contribution from the plumes is overall positive.Negative contributions originate from the environmental turbulence in the entrainment zone.Though the basic vertical structure of the flux contributions can be observed similarly over flat terrain (not shown), a noticeable difference is the area between x = 7 tkm and x = 9 km, where the moisture accumulations are reaching out from the ridges towards the valley center.Here, the turbulent moisture flux contributions are almost zero (Figure 3f,g).From this, it can be concluded that the horizontal moisture transport in this area is mainly achieved by the slope-wind circulation (see the horizontal component of the upper circulation cell in Figure 2b).

JPDFs over the idealized valley
Over the idealized valley, the JPDFs are calculated on horizontal levels of constant height (Figure 4).This means that different flow regimes are intersected.For instance, the 800 m level crosses the entrainment zone in the valley center and the slope-wind layer.At 400 m height in the valley center (Figure 4a,e), the JPDF still resembles the statistics over flat terrain.However, the range of the potential-temperature and specific-humidity fluctuations is clearly higher while the vertical-wind velocities are slightly smaller.At slope height (800 m, Figure 4b), the total turbulent heat flux is considerably diminished compared with the same height level over flat terrain.At 1,800 m (Figure 4d,h), both the turbulent heat and the turbulent moisture flux are nearly disappearing (i.e., the JPDFs are nearly centrally symmetric).

Turbulent heat flux
In order to represent the full vertical transport, the contributions from the slope-wind circulation can be included in the flux statistics (Figure 5).With that, the density function for the heat flux at 400 m (Figure 5a) is asymmetric towards positive temperature perturbations at low vertical wind velocities in contrast to Figure 4a, where a negative bulge was apparent.At 800 m (Figure 5b), the total heat flux only has a slightly smaller value compared than for flat terrain, although the influence of the coherent structures is considerably smaller.This leads to the conclusion that a large amount of heat is transported by the slope winds.Above, the JPDFs at higher altitudes are more difficult to interpret.This is because of the high horizontal inhomogeneity between the updraft regions over the ridges and the largely non-turbulent areas of subsidence over the valley center.Thus, the total vertical heat flux at these heights is generally small.

Turbulent moisture flux
The JPDFs of the vertical moisture flux including the contributions from the slope-wind circulation (Figure 6) show a similar picture.Primarily, the total flux values are significantly higher than over flat terrain and reach about twice the surface flux at 800 m (Figure 6b). Figure 6e-h indicates that a large amount of moisture is transported by the slope-wind circulation.Different to a flat surface, the , , 1, but for the reference experiment.The height level of (d) and (h) is raised to 1,800 m in order to cover the elevated regions over the ridges.cs: coherent structures.
, , coherent structures even include negative flux values in the valley (Figure 6a,b, fourth quadrant).This is explained by plumes of drier air advected from the valley center to the moist slope-wind layer by the horizontal mean flow.
The drying of the valley center is a result of the subsiding branch of the cross-valley circulation (Weinkaemmerer et al., 2022).A comparison of Figures 4h and 6d shows that the moisture transport by turbulence is nearly negligible at upper levels.

Role of the slope-wind circulation
Comparing Figures 5d,h and 6d,h shows that the thermally induced circulation also determines the shape of the total JPDFs as the turbulent perturbations are clustered around the local-circulation components.As already mentioned, the slope-wind circulation takes on a significant amount of the vertical heat and moisture transport.For the heat flux, this applies especially to , , F I G U R E 6 As in Figure 5, but for moisture.
the slope-wind layer (Figure 5f), whereas for the moisture flux the contribution from the mean circulation clearly exceeds the coherent-structure part, especially in the slope-wind layer (Figure 6f) and at the height of the hilltop plumes (Figure 6h).Generally, the maximum vertical mean wind values in Figures 5e-g and 6e-g can be assigned to the slope flow.Mostly, they also coincide with the highest  c and (q v ) c values.Conversely, data points with negative w c belong to the recirculation (see Figure 5g).In the case of the heat flux, this downward movement together with the upward movement of negative temperature fluctuations leads to a negative contribution to the heat flux (Figure 5g,h, second and fourth quadrants).For moisture, in contrast, the circulation part of the flux is dominated by the first and the third quadrants, meaning that its total value is overall positive.

Reference case
In this section, the simulations over three-dimensional terrain are analyzed.Figure 7a shows the result of the conditional sampling for the REF3D case at midday (1200 UTC) after 6 hr of simulation.The terrain features a group of hills with heights of 500 m above a surrounding plain.The sampled plumes are evenly distributed over the domain, and the orography does not seem to have a significant impact on their strength.The boundary layer is about 1.5 km deep at this time, meaning that the hills are rather shallow compared with the boundary-layer depth.In Figure 7a, artifacts of the subdomain decomposition are also visible, showing up as edges in the data field aligning with the subdomain boundaries.Again, the plumes cannot be treated independently from the thermal circulation.The effect of the slope-wind circulation becomes apparent in Figure 7c,e, showing that the temporally averaged plume height is maximal on the ridges.This impression is strengthened by the 12 hr average (Figure 7g).It means that the plumes are more frequent on the mountain tops, where they converge horizontally due to the upslope winds.This has a distinct effect on the time-averaged statistics and is in accordance with the positive vertical mean flow over the tops visible in Figure 3. Without time averaging, as in Figure 2c, no evident accumulation or intensification of the plumes over the hilltops was found.This is because the conditional sampling is based on the local tracer-concentration statistics (in particular, the standard deviation of the tracer concentration).Although seeming counter-intuitive at first glance, this is consistent with the decomposition of the flow into a domain average, a local mean value, and a turbulent fluctuation.With that, it can be seen that the scale of the mean circulation is sufficiently separated from the scale of the coherent structures.At the same time, this shows that the subdomains are reasonably sized.In other words, a low-elevation orography without steep gradients has no big impact on the instantaneous characteristics of the individual boundary-layer plumes as the turbulent motions live on a smaller scale than the thermal circulation.

Impact of an upper level wind
Similarly to the local mean flow, a large-scale, upper level wind can also be expected to have an impact on the temporally averaged plume density.In the WIND3D case, the direction of the prevailing background wind is northwesterly.As the comparison of Figure 7c,d shows, the effect of the large-scale wind only becomes apparent in the temporal average.The average over the first 6 hr of simulation shows that the frequency of the plumes on the windward, northwestern side of the ridges is diminished compared with REF3D.This is linked to an increase in the plume frequency on the leeward slopes.The differences intensify in the second half of the simulation (Figure 7e,f) when the plumes are generally less persistent on the wind-exposed hills compared with the surrounding plain.Owing to the generally taller plumes in the afternoon, this pattern also dominates the 12 hr mean (Figure 7g,h).
The upper level wind also influences the vertical structure of the boundary layer.The domain-averaged profiles of the horizontal wind components at 1800 UTC (and at 1200 UTC for comparison) are presented in Figure 8.The figure also shows the impact of the upper level wind on the potential-temperature and specific-humidity profiles compared with REF3D.Although the effect is small, the outcome confirms the finding of Weinkaemmerer et al. (2022) that the temperature profile stays almost unaffected whereas the export of moisture from a valley to the free atmosphere is reduced depending on the valley depth and the atmospheric stability.This also seems to be true over a three-dimensional orography of shallow hills without long, deep valleys.The 1200 UTC curves show that the upper level wind mainly takes effect in the second half of the day, at least in this case.This is when the mean distribution of the plumes differs most between REF3D and WIND3D (see Figure 7).At 1000 UTC, the domain-averaged flux profiles of REF3D and WIND3D still overlap (Figure 8e-g,i-k).This time has been picked to be representative for the late-morning conditions.Studying the contributions to the vertical moisture flux from turbulence and local circulation (Figure 8i-k) at 1600 UTC reveals that the moisture transport by the circulation is strongly reduced in the WIND3D case, which is partly counteracted by an increase in turbulent exchange, as described in Weinkaemmerer et al. (2022).This increase of the turbulent moisture flux traces back to the environmental turbulence and, partly, to the coherent structures (Figure 8l).Also, the in-structure variability plays a significant role in the WIND3D case.This indicates that the plumes are more irregular and inhomogeneous than in the reference case, where this contribution is almost negligible.The in-structure variability is also increased for the turbulent heat flux (Figure 8h).All in all, there is a strong sensitivity of the turbulent components and the local circulation to the upper level wind, especially for the moisture flux.This makes it difficult to draw a direct connection from the spatial plume prevalence (Figure 7) to the vertical moisture transport, as other turbulence and the slope-wind circulation generally play a big role.The main impact of an upper level wind on the vertical moisture transport is due to the reduced advective transport by the slope winds.

CONCLUSIONS
For this study, (semi-)idealized LESs have been carried out over both a quasi-two-dimensional periodic valley and a three-dimensional, smoothed real orography featuring a group of hills.Driven by either a prescribed, constant surface heat flux or by time-dependent solar radiation, an almost cloud-free CBL develops that is typical for daytime, fair-weather conditions.The CBL is characterized by thermal updrafts or plumes being part of the small-scale, turbulent motions.These plumes are often referred to as coherent structures or organized turbulence.In contrast to flat terrain, the boundary-layer structure over mountains and hills is horizontally inhomogeneous, and thermally induced slope winds emerge as a local mean flow.A conditional sampling method has been used to identify the coherent structures in order to analyze the simulations with a focus on the plume characteristics.The implementation of the conditional sampling is consistent with the flow decomposition presented in Weinkaemmerer et al. (2022), where a local average is defined besides a background value and a turbulent part.The main findings can be summarized as follows.
• It is found that, in convective conditions, the coherent structures in the mountain boundary layer are influenced by both the thermally induced winds and the background upper level wind.In the slope-wind layer, the plumes are moving upslope with about the speed of the slope wind.In consequence, the temporally averaged plume thickness is highest on the crests of the hills.The coherent structures make up the dominant contribution to the vertical heat flux.
• In analogy with the CBL over flat terrain, the non-local transport related to the boundary-layer plumes in the complex-terrain CBL is evident in the irregular shape of the JPDFs of vertical velocity, potential temperature, and specific humidity fluctuations.In the CBL over mountains, the features of the JPDFs can be well explained by the properties of the thermally driven mean circulation.The JPDFs show that the advection of coherent turbulent structures by the thermally driven flow is a dominant process in the CBL over orography.Generally, the coherent structures are less relevant here as the local circulation takes on a greater role.Also, the turbulent entrainment at the top of the CBL is partly replaced by the downward components of the mean circulation.
• An upper level wind reduces the vertical moisture transport and, especially on the windward hillsides, the frequency of plumes.The partitioning of the vertical moisture flux between plumes and local circulation is very sensitive to the upper level wind.Also, the non-organized turbulence outside of the plumes plays a more significant role if a large-scale wind is present.
For future research, the experimental set-up should be extended to various atmospheric conditions and different Bowen ratios.Also, the possibility of cloud formation should be taken into account in order to investigate how the vertical transport of heat and moisture by (organized) turbulence and mean flows and its sensitivity to a background wind changes with the onset of clouds.So far, the results of this study are also applicable to tracers or pollutants in general.
Regarding parametrization development, no clear path can be given.Several researchers suggest joint-Gaussian JPDFs to describe turbulent fluxes (Wyngaard and Moeng, 1992;Wang and Stevens, 2000;Chinita et al., 2018).However, the JPDFs over the idealized valley are far more asymmetric than over flat terrain.In fact, their shape is largely determined by the local circulation.This finding lends credibility to early attempts of parametrizing vertical heat transport over complex terrain on the basis of analytical models of the slope-wind system, such as the Prandtl model (e.g., Noppel and Fiedler, 2002).However, although the idea of parametrizing the mass-flux component of subgrid-scale fluxes with a slope-wind model is conceptually attractive, practical implementations must consider the documented shortcomings of the underlying theory.For instance, recent studies (e.g., Zardi and Serafin, 2015) demonstrated that the steady-state Prandtl model is heavily inaccurate at representing upslope wind speed and temperature anomalies at shallow slope angles and weak stratification.Also, the sensitivity of the individual contributions from turbulence and mean circulation to a background wind, especially for the tracer/moisture transport, complicates parametrization design.Even in high-resolution NWP simulations, where the slope winds are reasonably resolved, a parametrization should allow for the advection of subgrid turbulent structures, as we have shown that the advection of plumes is an important process in the CBL over mountainous terrain.
probability density functions of the turbulent vertical velocity fluctuation w ′ with (a)-(d) the turbulent potential temperature fluctuation  ′ and (e)-(h) the turbulent specific humidity fluctuation q ′v for a flat surface at 4 hr at different height levels (boundary-layer depth ≈1,500 m).Numbers in parentheses in the legends represent the total vertical flux and the contribution of the coherent structures (cs), each normalized by the surface flux.Iso-probability contours are at 0.05 × 10 −3 , 0.1 × 10 −3 , 0.2 × 10 −3 , 0.5 × 10 −3 , 0.001, 0.002, 0.005, and 0.01.Slope cross-sections for the reference experiment at 4 hr showing (a) the contours of the coherent structures (cs) obtained by conditional sampling and (b) the horizontal wind component of the slope-wind circulation and the isolines of the potential temperature (K).The mixed and the boundary-layer height are illustrated as a dashed and a solid line respectively.(c) Vertical thickness of the structures in top view and (d) horizontal component of the structures' bulk velocity at 10 m above ground.

E 3
Slope cross-sections for the reference experiment at 4 hr showing the (a)-(d) turbulent heat flux and (e)-(h) turbulent moisture flux.(a) Total turbulent heat flux; (b)-(d) the turbulent heat-flux contribution from (b) the coherent structures (cs), (c) the environment, and (d) the in-structure variability.(e) Total turbulent moisture flux; (f)-(h) turbulent moisture-flux contribution from (f) the coherent structures, (g) the environment, and (h) the in-structure variability.Also shown are the slope-wind vectors as well as (a)-(d) the potential temperature (K) and (e)-(h) the specific humidity (physical notation).

F I G U R E 5
As in Figure4a-d, but showing (a)-(d) the local-circulation components combined with the turbulent fluctuations and (e)-(h) only the local-circulation components of the total heat flux.cs: coherent structures.In (e)-(h), one data point stands for one intersected grid cell.Green shadings are close to the slope, and yellow ones are close to the valley center.
Vertical thickness of the coherent structures over three-dimensional orography (height of the isolines in meters) at 1200 UTC (after 6 hr simulation time) for (a) three-dimensional reference case (REF3D), (b) three-dimensional case with initial upper level wind (WIND3D); (c)-(h) the mean vertical thickness averaged over each subdomain at a 1 min interval over (c, d) the first 6 hr of simulation, (e, f) the second 6 hr of simulation, and (g, h) the whole simulation time.
Domain-averaged profiles of (a) the zonal wind u, (b) the meridional wind v, (c) the potential temperature , and (d) the specific humidity q v for the three-dimensional reference case (REF3D) and three-dimensional case with initial upper level wind (WIND3D) at initialization (0600 UTC) and different later times.(e) The total heat flux, (f) the turbulent heat flux, and (g) the local-circulation component.The turbulent part includes the subgrid contribution.(h) The decomposition of the resolved turbulent heat flux according to Equation (12) for REF3D and WIND3D; (i)-(l) analogously for the moisture flux.

TA B L E 1
List of simulation experiments.