Anabatic Flow along a Uniformly Heated Slope Studied through Large-Eddy Simulation

Abstract

Anabatic flows are common phenomena in the presence of sloping terrains, which significantly affect the dynamics and the exchange of mass and momentum in the low-atmosphere. Despite this, very few studies in the literature have tackled this topic. The present contribution addresses this gap by utilising high-resolved large-eddy simulations for investigating an anabatic flow in a simplified configuration, commonly used in laboratory experiments. The purpose is to analyse the complex thermo-fluid dynamics and the turbulent structures arising from the anabatic flow near the slope. In such a flow, three main dynamic layers are identified and reported: the conductive layer close to the surface, the convective layer where the most energetic motion develops, and the outer region, which is almost unperturbed. The analysis of instantaneous fields reveals the presence of thermal plumes, which are stable turbulent structures enhancing vertical transport and mixing of momentum and temperature. Such structures are generated by thermal instabilities in the conductive layer that trigger the rise of the plumes above them. Their evolution along the slope is described, identifying three regions responsible for the plumes generation, stabilisation, and merging. To the best of the authors’ knowledge, this is the first numerical experiment describing the along-slope behaviour of the thermal plumes in the convective layer.

1. Introduction

Anabatic and katabatic winds are frequent and widespread phenomena that characterise complex terrains [1,2,3,4,5,6,7]. They are triggered by the diurnal cycle, which continuously modifies the local temperature gradients and buoyancy along the topographic surfaces [1]. In the context of pure slope flows, katabatic winds are studied more due to their air quality implications. These winds typically arise at nighttime when the ground undergoes radiative cooling, generating a downslope flow and facilitating the return of pollutants or other substances downward, possibly into urban areas. Conversely, anabatic flows are buoyancy-driven flows initiated by solar heating of a slope [8] and sustained by the continuous (although heterogeneous) heat transport from the warm surface to the atmosphere [9]. They are more problematic to measure and simulate than katabatic winds, due to the higher vertical extension above the ground and the larger number of convective scales of motion that develops [10]. Despite being investigated less, the anabatic flows play a crucial role in the low atmosphere dynamics and impact the exchange processes between the lower and upper atmosphere. The onset of surface convection and development of the anabatic flows is crucial to break the flow stagnation at the slope base [11,12], thus favouring the pollutant dispersion. Moreover, transport and redistribution are enhanced through injection in larger mesoscale systems [13] or re-cycle to the slope base, as observed for heat transfer by Serafin and Zardi [14]. The heterogeneity in the anabatic flows modifies the structure of the convective boundary layer (CBL), causing transient depression in its depth [15] and more defined layered structures [16], and thus making the definition of the CBL height a matter of discussion [17].

The study of the anabatic flows has mostly relied on field or laboratory experiments to isolate the driving processes and improve our comprehension of their physics and has involved the use of analytical models based on the fundamental understanding of the developing mechanisms. An exhaustive review on hill flows is presented in Finnigan et al. [18]). One of the most accepted analytical models was proposed by Prandtl [19] to address the steady laminar flow over an infinite plate in a stratified atmosphere with Boussinesq approximation. Although rather simplified, the Prandtl model is used as a base for following the improvement of the theory, e.g., adding conditions of flow unsteadiness [20], heterogeneity [21], varying eddy viscosity [22], multiple dynamical balances [9] Coriolis force [23], periodicity [24] and weakly nonlinear effects [25]. On the other hand, pioneer field experiments have paved the way for key studies on the comprehension of the local processes driving the onset of the anabatic flows [26], the evolution of the anabatic layer [9] and its unsteadiness [27]. Laboratory experiments have delved into the initial development of the slope flow [11,28,29] and the implication of a later development on air quality [30]. Of more recent interest is the evaluation of rising plumes and flow separation from the slope due to the surface heating. The balance between the surface heating and the flow vorticity has been theorised by Crook and Tucker [31] at the base of the flow separation when the slope is uniformly heated, with local flow convergences (divergences) at the top (bottom) of the slope favouring vertical motions [14]. Further characterisation of the flow separation has been given by Princevac and Fernando [8] in terms of balance between surface heat flux, buoyancy force, and the plume size. Recently, laboratory experiments have performed similar investigations to address the slope-angle dependence of the flow separation. Hocut et al. [32] identified three regimes of separation according to the angle of a constantly inclined slope rising from a flat basin. For angles larger than 20${}^{\circ }$, the separation generates a plume fed by the upslope flow, while for angles smaller than 10${}^{\circ }$, the plume superimposes the slope flow. Within that angle range, the upslope flow requires the entrainment of an ambient flow to sustain the plume. Goldshmid et al. [33] revised these findings, adding a horizontal plateau at the top of the slope.

From the numerical perspective, few studies on this topic have been conducted using computational fluid dynamics (CFD) models and numerical weather prediction models. CFD simulations have mostly delved into the study of katabatic flows [34,35,36,37], whilst very little work has been carried out for anabatic flows. Schumann [38] performs a pioneering large-eddy simulation (LES) of an anabatic flow along a uniformly heated inclined plane, showing that the turbulence structures within the steady-state upslope boundary layer are a function of the slope angle and the surface roughness relative to a proper length scale of buoyancy. With a similar approach, Dörnbrack and Schumann [39] addressed the relation between orography and convective scales, while Anquetin et al. [40] investigate the anabatic-wind onset by analysing the formation and breakup of the cold pool over an idealised two-dimensional valley. Regarding the Reynolds-Averaged Navier–Stokes (RANS) approach, specific investigations of anabatic flows have not been reported in the literature, while their effects on the urban environment have been assessed by Luo and Li [41]. Direct Numerical Simulations (DNS) have been performed to reproduce idealised cases, such as uniformly heated infinite slopes or constant heat flow. Fedorovich and Shapiro [42] simulate katabatic and anabatic flows for different slope inclinations and temperatures, analysing the transition to turbulence and the possible parametrisation of mean profiles. Giometto et al. [43] perform several simulations changing the slope angle and temperature to derive information for improving current atmospheric models. In the framework of numerical weather prediction models, the Weather Research and Forecasting (WRF) model has been provided with a module to perform in LES-like mode [44]. Using this model, Catalano and Cenedese [45] were able to capture the diurnal cycle of slope winds and top-layer interaction with the valley breeze, while Serafin and Zardi [16] evaluated the diurnal heat transfer between slope and valley flows. Similar investigations were performed using other mesoscale models in LES-like mode giving a satisfactory slope-flow representation [46,47]. Furthermore, Weigel et al. [48] demonstrated the role of up-valley winds in stabilising the valley atmosphere despite strong thermal forcing at the ground.

The present work aims to gain knowledge about the basic mechanisms governing the anabatic flows. The focus is given to the turbulent structures, which eventually play a crucial role in the exchange processes at the surface. To this end, the anabatic flow developing in a simplified geometry is reproduced through a high-resolved LES. The configuration chosen reproduces a usual setup of laboratory experiments that investigate this topic: a uniformly heated double planar slope starting from neutral background stability. One of the added values of the numerical approach is the high space-time resolution of the dynamic fields, which enable overcoming the typical limitation imposed by the experimental measurement apparatus and gaining a more exhaustive comprehension of the specific process under investigation. To the best of the authors’ knowledge, this is the first high-resolved LES specifically designed to evaluate the turbulent dynamics of an anabatic flow in a simplified setup. From the numerical point of view, the reliable numerical simulation of complex thermo-dynamic interaction poses some challenges. First, thermal convective turbulence poses some stability issues and requires accurate numerical schemes to be resolved [49]. In addition, the transient nature of heat transfer is most effectively addressed by time-dependent simulations, which can reproduce the evolution of instantaneous configurations but require the use of a careful simulation approach and turbulence model.

The remainder of this paper is organised as follows: Section 2 presents the specifics of the LES model; Section 3 is devoted to the validation of the numerical approach; Section 4 introduces the case study of an anabatic flow over a uniformly-heated finite slope; Section 5 reports the results discussion. Finally, Section 6 draws the conclusions.

2. Numerical Methodology

The modelling approach for numerical simulations and the case settings are detailed in the following sections.

2.1. Governing Equations

The anabatic flows fall into the category of the thermally-driven flows, which are here modelled considering the fluid as incompressible and adopting the Boussinesq approximation to account for the onset of the buoyancy force. The thermal stratification of the fluid (if present) is described by the potential temperature, which is decomposed as $\theta ={\theta }^R+T$, where ${\theta }^R$ is the background temperature and T is its perturbation (see Prandtl [19]). Assuming that the background temperature is a linear function of the height z and $\beta$ is the thermal expansion coefficient, one can define the Brunt–Väisälä buoyancy frequency $N^2=\beta \left(d{\theta }^R/dz\right)$ as a constant. In this framework, the governing equations are written as follows:

$$\begin{array}{cc}\hfill \frac{\partial u_i}{\partial x_i}& =0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}& =-\frac{\partial p^{★}}{\partial x_i}+\nu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}-g_i\beta T,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\partial T}{\partial t}+u_i\frac{\partial T}{\partial x_i}& =-\frac{N^2}{\beta g}{u}_i{g}_i+\kappa \frac{{\partial }^{2}T}{\partial x_i\partial x_i},\hfill \end{array}$$

(3)

where $u_i$ are the velocity components, $p^{★}=\left(p-{\rho }_0{g}_i{x}_i\right)/{\rho }_0$ is the pressure deviation from the kinematic hydrostatic pressure, ${\rho }_0$ is the constant reference density, $g_i$ is the gravity acceleration vector and $g=9.81$ m/s${}^2$ is its modulus, $\nu$ is the kinematic viscosity, and $\kappa$ is the molecular thermal conductivity. It is worth noticing that Equations (2) and (3) simplify if expressed in terms of buoyancy force modulus $b=g\beta T$ (where T is the usual difference of temperature), which is convenient in specific applications [43,50]. In absence of thermal stratification ${\theta }^R=constant$ and $N^2=0$.

2.2. Non-Dimensional Numbers and Parameters

Two non-dimensional numbers have been selected to describe the natural convection in the presence and absence of thermal stratification. Following Giometto et al. [43], the Grashof number $Gr$ for the thermally stratified flows is defined as:

$$Gr=\frac{|b_S{|}^4}{{\nu }^2{N}^6}=\frac{g^4{\beta }^4{T}_S^4}{{\nu }^2{N}^6}$$

(4)

where $T_S$ is the temperature perturbation at the slope. This represents the ratio of the buoyancy $b_s=g\beta T_S$ to the viscous and stratification forces. In absence of thermal stratification, the Rayleigh number

$$Ra=\frac{g\beta |T_s|\mathcal{H}}{\nu \kappa }$$

(5)

describes the regimes of natural convection, being the ratio between momentum diffusivity and thermal diffusivity. Here, $\mathcal{H}$ is the characteristic length of the system. Additional parameters are the molecular Prandtl number $Pr=\nu /\kappa$ (ratio of momentum diffusivity to thermal diffusivity) and the buoyant velocity $U_b=\sqrt{g\beta |T_s|\mathcal{H}}$, which is the convective velocity rising at a heated/cooled surface.

2.3. Large-Eddy Simulation Technique

Simulations are carried out by adopting an LES approach: the scales of motion larger than the computational grid are directly resolved, while the smallest are modelled through a sub-grid scale (SGS) model. Among the others, see Sagaut [51] for an overview of this topic. As known, small scales exhibit more universal features compared to the large scales: they are more homogeneous, isotropic, low in energy and only responsible for dissipating turbulence kinetic energy. Hence, the modelling of such scales can be more effective with respect to the turbulent models that have to account for the entire range of scales of turbulence, like those used for (Unsteady) RANS simulations (see Pope [52]).

Numerically, the small scales are implicitly determined by the computational grid. The latter acts to the velocity field as an implicit spatial filter (denoted with an overlying bar), which leads to a decomposition into a filtered and a residual component: $u_i=\overline{u_i}+u_i^{\prime }$. The former represents the resolved large scales, the latter accounts for the sub-grid scales to be modelled. Analogous decomposition applies to the other fields. Filtering the momentum Equation (2), the SGS stress tensor ${\tau }_{ij}$ appears

$$\frac{\partial {\tau }_{ij}}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\overline{u_i{u}_j}-{\overline{u}}_i{\overline{u}}_j\right),$$

(6)

and needs to be modelled. We adopt the Smagorisky model [53] to close the system, where the isotropic part of the SGS tensor ${\tau }_{kk}{\delta }_{ij}/3$ is included in the pressure term, while the deviatoric part reads:

$${\tau }_{ij}-\frac{{\delta }_{ij}}{3}{\tau }_{kk}=-2{\nu }_{sgs}\overline{S_{ij}}=-2c_s^2{\overline{\Delta }}^2\left|\overline{S_{ij}}\right|\overline{S_{ij}},$$

(7)

where the sub-grid scale viscosity ${\nu }_{SGS}$ is implicitly defined, $c_s^2=0.4225$ is the Smagorinsky constant, $\overline{\Delta }={\left({\Delta }_x{\Delta }_y{\Delta }_z\right)}^{1/3}$ in the estimation cell width and ${\Delta }_i$ is the extension of the cell in the i-direction, and $\overline{S_{ij}}=(1/2)[(\partial {\overline{u}}_i/\partial x_j)+(\partial {\overline{u}}_j/\partial x_i)]$ is the resolved strain rate tensor and $|\overline{S_{ij}}|={\left(2\overline{S_{ij}}\overline{S_{ij}}\right)}^{1/2}$ is its magnitude. The Van Driest function [54] is also applied near the solid boundaries to correctly damp the SGS viscosity within the wall boundary layer.

Similarly, filtering of the temperature Equation (3) generates an SGS thermal flux ${\varphi }_i$, which is modelled by the flux-gradient hypotheses using:

$${\varphi }_i=\left({\overline{Tu}}_i-\overline{T}{\overline{u}}_i\right)=-{\kappa }_{sgs}\frac{\partial \overline{T}}{\partial x_i};$$

(8)

and the sub-grid scale thermal diffusivity is expressed by ${\kappa }_{sgs}={\nu }_{sgs}/P{r}_t$, i.e., by applying the Reynolds analogy assumption and setting the turbulent Prandt number to $P{r}_t=0.7$. This assumption has been proven to be effective in the simulation of natural convection processes [55,56].

2.4. Algorithm and Numerical Settings

The software used for the simulations is OpenFOAM [57], which is an open-source toolbox written in C++. Several models and algorithms are already implemented, but the resolution of Equations (1)–(3) requires an ad hoc numerical code that has been developed on the basis of the existing pimpleFoam solver. The solver is based on the Pressure-Implicit with Splitting of Operators (PISO) algorithm.

Overall, second-order schemes are used for time and space discretisation. Time derivative is computed with an implicit Euler backward scheme, whilst the central-difference scheme is used for all spatial derivatives except the advective terms in the temperature equation. This is discretised by the Gauss–Gamma scheme (${\gamma }^{*}=0.2$) to stabilise the computation. Such scheme is a bounded central-difference scheme based on the normalised variable diagram [58]. The time step is dynamically computed to enforce the Courant–Friedrichs–Lewy condition $CFL<0.5$, which ensures the stability and accuracy of the simulation [59]. A numerical correction is applied to prevent non-physical values of T. Such values can be triggered by numerical oscillations arising in the early stages of flow development and involve a limited number of cells that tends practically to zero once the statistical steady-state is reached.

3. Simulation Approach Validation

The numerical solver and settings are validated by reproducing the paradigmatic case of an infinite heated slope that is sketched in Figure 1. It consists of a heated and infinitely long plane, inclined at an angle $\alpha$ to the horizontal (see Figure 1b). For the sake of simplicity, the flow is described in the system of reference of the slope: $\zeta$ and $\xi$ are the slope-normal and slope-parallel direction, respectively. The anabatic motion is driven by the buoyancy force generated by the temperature gradient at the slope surface, which has two components along and normal to the slope. Such a configuration is based on the milestones of Prandtl [19]. It allows for a precise analysis of the fundamental physical phenomena underlying anabatic flows, generated by the nonlinear thermo-fluid dynamic interaction, without the additional complications given by complex geometries.

3.1. Validation Case Setup

The validation is made against the DNS by Giometto et al. [43], which imposes a constant temperature at the slope and a background thermal stratification. Particularly, the case U30H is reproduced and the results are expressed in terms of buoyancy field b rather than temperature variation. The flow regime corresponds to $Gr=2\times {10}^{11}$ with a slope inclination of $\alpha ={30}^{\circ }$ and the air viscosity $\nu =1.5\times {10}^{-5}$. Quantities are made non-dimensional by means of system characteristic length $L=b_S/N^2$, velocity $U=b_S/N$, and buoyancy $B=b_S$. The non-dimensional time is $\mathcal{T}=2\pi /sin\alpha$, which defines a characteristic time period.

The computational domain is shown in Figure 1a and is three-dimensional, including the (spanwise) y-direction. Its extension in the non-dimensional unit is $0.241\times 0.241\times 0.324$ in the $\xi$-, y-, $\zeta$-direction, respectively. The mesh counts $96\times 96\times 256$ cells, which are stretched near the slope surface by a hyperbolic tangent function [60] to ensure that the wall-boundary layer is directly resolved. To this end, the first cell near the slope has the same height as the cell in the reference DNS.

To simulate an infinite slope, cyclic conditions are imposed at the boundaries perpendicular to the slope. The top boundary simulates a free interface (zero-gradient condition for all the variables). At the slope, a non-slip condition for velocity $u_s=0$ m/s, constant positive surface buoyancy $b_s=1$ m/s${}^2$ and zero-gradient for other variables is imposed. A sponge region of zero gravity force $g_i=0$ m/s${}^2$ is imposed in a region corresponding to the $15\%$ of the top of the domain to avoid possible numerical instability, as suggested in Giometto et al. [43].

The simulation is run for $6\mathcal{T}$ time periods and the statistics are collected in the last $3\mathcal{T}$, when the statistical steady-state configuration is reached. Quantities are averaged along the streamwise ($\xi$) and spanwise (y) directions, as well as in time. The root-mean-square (RMS) are computed as ${\psi }_{rms}={\langle {\left[\psi -\langle \psi \rangle \right]}^2\rangle }^{1/2}$, where $\psi$ is a generic variable and the angular brackets denote the space-time average.

3.2. Validation Results

First- and second-order statistics from LES are reported and compared with the reference DNS [43]. The analyses are performed in the system of reference of the slope; thus, the profiles are taken along the slope-normal direction.

Figure 2a shows the non-dimensional averaged profiles of streamwise velocity and buoyancy, as well as the RMS of velocity components and buoyancy magnitude. The LES mean profiles practically collapse onto the DNS profiles. The velocity RMS components are slightly underestimated by the LES in a region very close to the slope, where the convective motion is more energetic. The buoyancy RMS reproduces the peak at $\zeta /L={10}^{-3}$ even if its intensity is slightly overestimated.

Figure 2b displays the slope-normal profiles of non-dimensional turbulent kinematic momentum and buoyancy fluxes, and the resolved turbulence kinetic energy $k=1/2\left(u_i^{\prime }{u}_i^{\prime }\right)$. Turbulent fluxes are in very good agreement with the reference profiles, while the LES turbulence kinetic energy exhibits a limited underestimation consistent with RMS profiles.

Overall, the accuracy of the mean and turbulent fields is satisfactory, despite minor discrepancies in the RMS profiles. Thus, the present simulation approach allows to effectively reproduce the thermo-fluid dynamic of the system. The LES results are then applied for simulating anabatic flows in the case study described and analysed in the following sections.

4. Case Study Setup

The anabatic flow along a finite slope uniformly heated is studied in a typical geometry of a laboratory experiment [28,32,61]. Figure 3b sketches the case setup: a double slope is placed in the centre of a water tank with dimensions of $2.0\times 0.6\times 1.2$ m, respectively, for width (x), depth (y) and height (z). The slope is $L=0.3$ m long, inclined of an angle $\alpha ={35}^{\circ }$ with respect to the ground, and heated up of $T_s=2.78$ K compared to water. Such difference of temperature has been chosen to generate an average heat flux at the slope equal to $Q=q/A_s=900$ W/m${}^2$, where $q=\kappa C_p{\rho }_{0}dT/dn$ is the heat flux density at the slope, $dT/dn$ the temperature normal gradient, $C_p$ is the water specific heat at constant pressure, and $A_s$ is the area of the slope. The physical parameters for water are: thermal expansion $\beta =2.1\times {10}^4$ K${}^{-1}$, laminar and turbulent Prandtl numbers $Pr=6.13$ and $P{r}_t=0.7$, respectively. The boundaries of the tank are solid walls, except for the front ($y/H=0$) and back ($y/H=2$) vertical boundaries to which a cyclic boundary condition is imposed to make the domain virtually infinite along the direction of homogeneity (y).

It is worth highlighting that the configuration chosen for the present work is representative of experimental studies in a water tank, which are characterised by low $Ra$ compared to the real-system configuration. Schumann [38] showed that atmospheric flow can be reproduced by using wall functions, i.e., modelling the near-ground physics to reduce the computational cost of the simulation. This allows the simulation of high-$Ra$ flows but does not directly resolve the slope-boundary layer and the underlying physical mechanisms, which are the focus of the present work. For this reason, the study is concerned with highly resolved flow at low $Ra$.

The domain height and width are large enough to ensure that thermal convection along the slope reaches a configuration of equilibrium for a period of time sufficiently long to collect statistics. During this period, the influence of the boundaries is limited to the upper part of the tank. The analyses are focused on an observation region surrounding the slope, not significantly affected by the overall tank circulation (see discussion of Figure 4). The extension of the observation region is $1.5L\times 1.5L$ in the x- and z-directions and the entire domain in the y-direction.

The characteristic length is chosen to be the height of the slope $H=Lsin(\alpha )$ and the characteristic velocity is the buoyant velocity $U_b$. In such a way, the characteristic velocity is higher when the slope is more inclined and reduces to zero when it is horizontal (i.e., when no anabatic flow develops). Consequently, the characteristic time scale is $\mathcal{T}=H/U_b$ and the Rayleigh number of the system is $Ra=3\times {10}^8$, indicating a fully turbulent regime.

4.1. Computational Grid, Initial and Boundary Conditions

Taking advantage of the symmetry of the problem along the vertical plane at the mid-width (see Figure 3b), the computation domain is reduced to half of the tank and a symmetry boundary condition is applied to virtually simulate the remaining part. This reduces the computational cost without losing information about the physics of the problem.

Figure 3a displays the numerical mesh, which is structured and staggered, composed of $192\times 96\times 192$ cells in the x-,y-,z-direction respectively. A hyperbolic tangent stretching [60] is applied to vertical (z) and horizontal (x) directions in order to increase the resolution near the slope and the symmetry plane, where the more energetic motion is localised. In the vertical direction, the smallest cell at the slope surface is $6\times {10}^{-4}$ m high. Preliminary simulations have been run to state that the slope-boundary layer is directly resolved, i.e., ensuring that the first computational point at the slope is within $z^{+}=u_{\tau }z/\nu <1$ (where $u_{\tau }$ is the friction velocity) [62]. Such preliminary simulations have been performed on several meshes with an increased resolution above the slope and increasing stretching at the slope surface. Particularly, three meshes have been tested with an overall increase of about $20\%$ of cells in all directions above the slope (in the observation region). The present mesh is the intermediate between the three since the additional increment of points does not change the first- and second-order statistics.

In order to discuss the grid quality, as well as to assess the range of turbulence scales captured by the computational grid, the power spectra are reported [63,64]. Figure 5 shows the power spectra of velocity components and temperature perturbation in the system of reference of the slope. Data are collected in time at three selected points, which are placed at the mid-slope and at increasing distance from the surface (within the region where the most energetic anabatic flow develops, see also Figure 6). Velocity power spectra exhibit the $-5/3$ slope characteristic for the inertial range, while the steeper slope is associated with the inertial-diffusive range for buoyancy-driven flows [63]. In the temperature perturbation spectra, a decay of slope $-3$ is detectable. These verified that the inertial range is covered by the computational mesh adopted, which is a requirement for performing an LES. In addition, the value of ${\nu }_{SGS}$ is compared with the molecular kinematic viscosity $\nu$, by scrutinising the ratio $\gamma =\nu /(\nu +\langle {\nu }_{SGS}\rangle )$. Such parameter measures the contribution of the sub-grid scale with respect to the computational mesh used: $\gamma \sim 1$ indicates that the LES reduces to a quasi DNS, $\gamma \sim 0$ means a high contribution of turbulent model and thus a coarse spatial resolution, while $\gamma >1/3$ corresponds to $\nu /{\nu }_{SGS}=2$ and suggests a good mesh resolution [55,64]. In the region above the slope where the largest part of the anabatic flow is localised, $\gamma >0.75$ indicating that the correction due to the turbulent model is limited and that the turbulent kinetic energy is dissipated mainly by the resolved scale of motion than the sub-grid scales. Hence, the spatial resolution can be considered satisfactory.

The boundary conditions are: at the ground and at the slope, the velocity enforces the non-slip condition, pressure and turbulent quantities have zero-gradient. The temperature is set to constant value $T_s$ at the slope and to reference temperature ${\theta }^R$ elsewhere. At the top and lateral walls, the Spalding wall function [65] is applied to velocity and pressure, while turbulent quantities enforce the zero-gradient condition and temperature is constant ${\theta }^R$. Wall functions are applied since the analyses are not focused on the far-slope dynamics and high accuracy is not required. As mentioned above, the symmetry condition is applied to all variables at the symmetry plane [57]. The front and back boundaries are cyclic. As the initial condition, the water is neutrally stratified ($N=0$), at the constant temperature ${\theta }^R=296.15$ K, and at rest ($u_i=0$).

4.2. Overall Circulation

The simulation is run for more than $180\mathcal{T}$. The velocity time-series at selected locations in the proximity of the slope (not reported) indicate that the flow reaches a stable configuration after $70\mathcal{T}$ from the initialisation. The statistics are collected for a period of $110\mathcal{T}$ after reaching the stable configuration. Quantities are averaged in time and space, along the spanwise direction.

The overall circulation within the tank is briefly described. Figure 4a shows the non-dimensional mean velocity magnitude and streamlines across the entire computational domain. The circulation is driven by natural convection generated at the heated slope. The flow increases its intensity from the slope base to the apex, where it separates from the surface and gives rise to an ascending flow. A principal vortex dominates the top-tank zone ($z/H>2.0$), while the bottom-tank zone ($z/H<2.0$) is not impacted by the principal vortex. At the vortex centre, the fluid is almost at rest. The hot fluid is carried downward from the top-tank zone by the principal vortex and comes in contact with the cold fluid in the bottom-tank zone (on the opposite side to the slope location). The average temperature of the top-tank is increased by the injection of hot fluid, transported by the ascending flow and mixed by the principal vortex. Instead, the bottom part practically maintains the initial temperature (see Figure 4b). The result is that the more energetic part of the descending flow from the top-tank deviates horizontally in the range $2.0<z/H<3.5$. Above the slope, the principal vortex bends towards the slope. However, the weak velocity intensity of this downward flow does not substantially influence the highly energetic anabatic flow near the slope (see also the discussion of Figure 6 and Figure 7). In the following Section 5, the focus is restricted to the observation region to detail the thermal and dynamic characteristics of the anabatic flow.

5. Results and Discussion

The anabatic flow along the slope is investigated by analysing the first- and second-order statistics distribution at a vertical plane normal to the spanwise direction (y) and along selected lines perpendicular to the slope surface. Three-dimensional and instantaneous turbulent features are also visualised and described.

5.1. First-Order Statistics

Figure 7 shows the spatial distribution of non-dimensional mean velocity components and temperature in the observation region. The higher values of horizontal velocity $u_x$ are mainly localised in a limited layer above the slope, where the anabatic flow develops. Hereafter, we refer to this region as the convective layer. The thickness of such a layer is almost zero at the slope base and increases approaching the apex, before being deteriorated by the ascending flow. It does not exceed $0.6H$ in the slope-normal direction, where the velocity and temperature perturbation profiles at $\zeta =2/6,3/6,4/6$ (that are less affected by the flow conditions at the of the slope edges) practically collapse on the same profile. An accurate definition of the layer thickness development is not reported since it requires a dedicated analysis that is beyond the aims of the present study. The streamlines show that the anabatic flow produces a weak horizontal flow that drives the fluid close to the ground towards the slope. The vertical velocity reaches the highest values above the slope apex ($z/H>1.2$), which is dominated by the ascending flow. Low positive values are also detectable in the convective layer. The streamlines highlight the presence of a weak downward flow from the top-tank zone (as discussed in Section 4.2). Above the centre of the slope, the flow is almost slope-normal directed, while being deviated in slope-parallel direction by the anabatic flow. The temperature distribution shows that the larger temperature perturbation is located extremely close to the slope surface. The contour lines follow the mean velocity field, except for the line at $T/T_s=0.01$ that deviates horizontally and stabilises at height $z/H=2.0$ in the bottom-tank zone. This confirms that the fluid is not stratified at the bottom part of the tank, while an unstable stratification increases climbing the slope.

Figure 6 shows the profiles of the non-dimensional mean velocity components (slope-parallel $u_{\xi }$ and slope-normal $u_{\zeta }$) and temperature, taken at five slope-normal lines equally distributed from the base to the apex of the slope. The profiles of the slope-parallel velocity $u_{\xi }$ are characterised by a rapid increase of velocity near the slope until a peak is reached, and subsequent decay to an edge-point above which the velocity shows a quasi-linear behaviour. As the apex of the slope is approached, the velocity exhibits a larger peak, which is displaced further away from the slope. Consistently, the region characterised by positive velocity and the edge-point increases in height. The profile at level $\xi /L=5/6$ shows the influence of the ascending flow, which prevents the growth of the velocity peak and transport momentum towards the slope-normal direction. The slope-normal velocity $u_{\zeta }$ profiles show negative values in the convective layer ($\zeta /H<0.5$), except at level $\xi /L=5/6$ where the ascending flow dominates. These velocities are one order of magnitude smaller than the slope-parallel ones. Negative velocities are due to the presence of the ground, and they are not present in the infinite slope case (see Figure 2a): the flow that reaches the slope from the left side of the bottom-tank is directed parallel to the ground; thus, the velocity that impacts the slope has a negative component in the slope-normal direction. The local minimum of velocity has an almost constant value $u_{\zeta }^{min}\sim 0.1U_b$, and its elevation approximately coincides with the height of the edge-point of slope-parallel velocity. In the surrounding area ($\zeta /H>1.0$), the vertical velocity is expected to tend to zero, but negative values are detectable. They are mainly generated by the downward flow from the tanks upper part. Temperature perturbation displays an exponential-decay layer ($\zeta /H<0.01$), where all the profiles decay exponentially below $T/T_s<0.15$ and practically collapse to the same curve (estimation done on non-reported profiles in logarithmic scale). This is a conductive layer where conduction dominates over convection. Above, there is the convective layer, where temperature experiences an almost linear decrease. The temperature maintains higher values as the apex is approached. Far from the slope surface ($\zeta /H>0.5$), the temperature field is almost unperturbed.

5.2. Second-Order Statistics

Figure 8 presents the distribution of turbulent kinematic momentum and heat fluxes at a vertical plane in the observation region. The momentum flux $\langle u^{\prime }{w}^{\prime }\rangle$ exhibits large positive values close to the slope surface and also in the final part of the convective layer $x/H>5.1$ (corresponding to $\xi /L>0.5$), with a peak in the vertical range $1<z/H<1.2$. This is the area where the downward flow from the upper part of the tank came in contact with the upper boundary of the convective layer and is deviated horizontally (see Figure 7), producing a positive spot of turbulent flux. The vertical heat flux $\langle w^{\prime }{T}^{\prime }\rangle$ assumes positive values in the convective layer, with a climax at the slope surface and above the apex where the ascending flow rises. They denote the tendency of turbulence to diffuse temperature upward within the convective layer.

Figure 9 displays the profiles of turbulent kinematic momentum and heat fluxes, as well as the turbulence kinetic energy $k=0.5u_i{u}_i$ and temperature perturbation RMS along selected slope-normal lines, in the system of reference of the slope. The turbulent momentum flux $\langle u_{\xi }^{\prime }{u}_{\zeta }^{\prime }\rangle$ has a maximum that corresponds to the inflexion point of the mean velocity component profiles in Figure 6, as predicted by the gradient-diffusion hypothesis. Above the maximum, all profiles rapidly decrease and reach an almost constant value approximately at the same height of the edge-point of the slope-parallel velocity. Hence, positive fluxes are confined in the convective layer. As expected, the intensity increases as the apex is approached and the main flow became more energetic. Turbulent slope-parallel heat fluxes $\langle u_{\xi }^{\prime }{T}^{\prime }\rangle$ show a strong peak within the conductive layer, where all the profiles practically overlap. In the convective layer, fluctuations have an almost linear decrease and assume zero values in the outer region. Turbulent slope-normal heat fluxes $\langle u_{\zeta }^{\prime }{T}^{\prime }\rangle$ peak around $\zeta /H=0.1$, where the mean temperature profiles have inflexion points. They also tend to zero in the outer region. Both sets of vertical and horizontal profiles indicate that turbulent diffusion is limited in the convective layer and acts to diffuse temperature in slope-parallel and -normal directions. The positive values of turbulence kinetic energy $\langle k\rangle$ are confined near the slope. Moving from the base to the apex of the slope, the intensity of k and the elevation of the profile peak increase. This is in accordance with the turbulent momentum flux and the slope-parallel velocity component. The profiles of temperature RMS $\langle T_{rms}\rangle$ exhibit a maximum at $\zeta /H=0.01$, i.e., the limit of the conductive layer. This is expected because of the strong thermal gradient arising at the interface. Above, these profiles rapidly decrease to zero and collapse almost onto the same curve within the convective layer, except for the line $\zeta /H=1/6$ that has an even more rapid decay.

Overall, three main regions of motion are detectable and qualitatively described in terms of dynamic and thermal characteristics:

Conductive layer: 

The thin region close to the slope surface ($\zeta /H<0.01$) where the slope-parallel velocity rapidly increases, temperatures following an exponential decay, and slope-parallel heat fluxes turbulent flows reach their maximum. Slope-normal heat transfer is dominated by conduction.

Convective layer: 

The layer above the slope where the anabatic flow mainly develops, corresponding to the region of largest slope-parallel velocity component and most energetic turbulent fluxes. Temperature profiles display an almost linear decrease, and the layer thickness increases along the slope before being destroyed by the ascending flow.

Outer region: 

Relatively far from the slope ($\zeta /H>0.6$). The velocity loses most of the energy generated by the buoyancy force at the slope, the temperature perturbation is almost negligible, and turbulent fluxes decay to low values. This region is almost unperturbed by the slope anabatic flow.

As already mentioned, the convective layer is still in evolution and does not reach a stable thickness in this configuration. For better analysis and quantification of such a layer, a more extended slope or a case study at higher $Ra$ would be required.

5.3. Turbulent Structures

Taking advantage of the opportunity given by the LES to reproduce instantaneous turbulent fields, the turbulent structures arising in the anabatic flow are analysed here.

Figure 10a shows the non-dimensional instantaneous temperature perturbation at the plane $y/H=0.3$, i.e., vertical and perpendicular to the slope. The temperature distribution reveals the presence of turbulent thermal plumes within the convective layer. Such plumes are generated near the slope base by thermal instabilities arising in the proximity of the surface. Subsequently, they grow vertically as they are transported along the slope by the anabatic flow. Near the slope apex, the plumes start to deteriorate and are finally destroyed by the ascending flow.

Figure 10b displays the three-dimensional structures of the turbulent thermal plumes through the visualisation of the isothermal surface $T/T_s=0.1$. Plumes emerge near the slope base as a sort of high-temperature bubble, which rapidly increases their height. At the centre of the slope, they are elongated in the vertical direction and presents a cylindrical base that is curved in the slope-parallel direction at the top. This hook-like shape is due to the main flow, which transports the plumes parallel to the slope. Near the apex of the slope, it is worth noticing that the axis of the plumes rotates to be almost horizontal before being destroyed by the ascending flow.

Figure 10c reports the temperature distribution near the slope surface (distance $\zeta /H=0.02$) to better understand the mechanism that triggers the rise of a thermal plume. Spots of high temperature are clearly detectable in the proximity of the slope surface, at the boundary between the conductive layer and the convective layer. They can be interpreted as a periodic perturbation of the conductive layer, and a careful visualisation of the plumes reveals that there is a thermal spot at the base of each of them. Such thermal spots are possibly generated when the cold flow from the bottom-tank impacts the hot slope surface. The strong unstable thermal gradient produces thermal instabilities (similar to the Taylor instabilities) in the conductive layer. Such instabilities give rise to the thermal plumes and are transported along the slope. Near the apex, the plumes detached from the thermal spots: the former are rotated horizontally by the main flow before being destroyed by the ascending flow; the latter merge together and produce structure elongated in the spanwise direction in the proximity of the slope apex.

Figure 11 summarises the structure and the evolution of turbulent thermal plumes. Their evolution along the slope can be divided into three main regions:

• Generation region ($0<\zeta /L<0.2$): Spots are almost absent at the bottom edge of the slope. The weak thermal instabilities are generated, but plumes are almost not detectable.
• Stabilisation region ($0.2<\zeta /L<0.6$): Spots assume an almost circular shape and the central temperature increases. Plumes emerge and develop in the convective layer, assuming their hook-line shape.
• Merging region ($0.6<\zeta /L<1.0$): Spots merge together while they are transported near the apex, where they appear elongated spanwise and are destroyed. Plumes detached from the surface and their axis is rotated horizontally. They are transported by the main flow and eventually destroyed by the highly energetic ascending flow.

Thermal plumes are turbulent structures confined within the convective layer and mainly develops vertically, thus transporting temperature from the slope surface across the convective layer. They can be responsible for the high level of turbulent thermal fluxes (see Figure 8) and contribute to the temperature mixing, which results in the almost linear profile of the temperature in the convective layer. Therefore, such structures play a crucial role in the dynamics of the convective layer and in the processes of vertical transport and mixing.

6. Conclusions

The high-resolved large-eddy simulation (LES) methodology is used to study an anabatic flow in a simplified geometry consisting in a uniformly heated double slope in an isothermal water tank. The configuration chosen reproduces the common setup of laboratory experiments. The Rayleigh number of the flow is $Ra=3\times {10}^8$ (fully turbulent regime) and the simulation is three-dimensional and time-dependent, enabling the analysis of instantaneous turbulent characteristics. To the best of the authors’ knowledge, this is the first LES specifically designed to evaluate the turbulent dynamics of an anabatic flow in a simplified setup.

The numerical approach is successfully validated against direct numerical simulations of an idealised case. Subsequently, the first- and second-order statistics of the selected case study are discussed. The overall tank dynamics shows that the anabatic flow generates a principal recirculation vortex, which is confined in the top part of the tank and does not significantly influence the development of the anabatic flow. The dynamics in the proximity of the slope are characterised by the presence of three characteristic regions: (i) the conductive layer, very close to the hot surface where temperature decreases exponentially and velocity rapidly increases; (ii) the convective layer dominated by convection, where the temperature has an almost linear decrease with the distance from the slope, and the most energetic mean and turbulent motions are confined; (iii) the outer region, which is almost unperturbed.

The analysis of instantaneous fields reveals the presence of thermal plumes, sketched in Figure 11. They are periodic turbulent structures continuously generated in the proximity of the slope base, transported along the slope and destroyed near the apex. These plumes arise from the vertical motions triggered by the thermal instabilities at the boundary between the conductive layer and the convective layer. Thermal plumes are confined into the convective layer and enhance vertical transport of momentum and heat, increasing the turbulence intensity and the mixing within the layer.

The present contribution gives the first description of such structures, which are described with respect to the anabatic flow dynamics. Additional analyses will be reported in future works.