Interactive comment on “ Adaptation of the meteorological model Meso-NH to laboratory experiments : implementations and validation ” by

Abstract. A meteorological numerical model can be a powerful tool to complement laboratory experiments applied to atmospheric sciences, but it needs to be adapted in order to represent flows generated in laboratory. This paper presents such an adaptation of the atmospheric non-hydrostatic model Meso-NH. The usually neglected viscous diffusive fluxes have been added to the model’s equations, along with the coding of an explicit bottom no-slip boundary condition was implemented. These implementations are validated against exact solutions of ideal flows. Meso-NH is then used in a configuration matching a laboratory experiment performed in the CNRM large stratified water flume meant to reproduce a neutral atmospheric boundary layer. The model is run for the first time in an explicit mode (Direct Numerical Simulation – DNS) at a very high resolution (1 mm) over a large grid. The comparison with the experimental data shows that the boundary layer height and vertical profiles of mean velocity are well captured by the model. This result further validates the implementations carried out in Meso-NH which can now be used in a DNS mode to simulate channel flows. The joint use of Meso-NH and laboratory experiments, along with the possibility to run DNS with Meso-NH, could lead to new findings or improvements in the field of atmospheric sciences.


Meso-NH
Meso-NH (Lafore et al., 1998) is a non-hydrostatic mesoscale research model developed by the Laboratoire d'Aérologie (UMR5560 UPS & CNRS, France) and the CNRM (UMR3589 Meteo-France & CNRS). It simulates atmospheric phenomena ranging from the synoptic scale to the large turbulent eddy scale, and it can be applied to real situations as well as aca-20 demic cases. Meso-NH is a grid point model using the Arakawa-C grid and the Gal-Chen and Somerville (1975) vertical coordinate ( Figure 1). The model's dynamics is based on the anelastic approximation (filtering of the acoustic waves). The prognostic variables are the three components of the wind, the potential temperature, the mixing ratio of hydrometeors and the turbulent kinetic energy -because of the anelastic approximation, the pressure is a diagnostic variable. The numerical schemes are an eulerian fourth-order centred advection scheme or a Weighted Essentially Non Oscillatory (WENO) scheme 25 of fifth-order (Shu and Osher, 1988) for the momentum components associated respectively with a leapfrog or a third order Runge-Kutta time marching (Lunet et al., 2017) and the Piecewise Parabolic Method (PPM) (Collela and Woodward, 1984) advection scheme for scalar variables like temperature, associated with a forward-in-time marching. Three types of lateral boundary conditions are available: rigid wall, cyclic and open boundaries. Meso-NH can be coupled to the SURFEX surface model (Masson et al., 2013) and it includes a large set of parametrizations. The model is freely available and can be down-30 loaded at http://mesonh.aero.obs-mip.fr.
In the present study, we use the v5-3-0 version, and the configuration designed for academic cases with dry air. We work under the Boussinesq's hypothesis in order to take into account the incompressibility of water. This is necessary since our main objective is to reproduce laboratory experiments using water as a fluid. Thus, we have a single phase incompressible fluid, as it is the case in a hydraulic flume. The Earth's sphericity is neglected and the horizontal coordinate system is cartesian. We also neglect the Coriolis force which does not affect the flow at the scale of a laboratory. The coupling with the surface model is disabled since we simulate laboratory experiments. No paramatrizations are activated and the implicit diffusion scheme is also turned off. In other words, Meso-NH is run explicitly.

5
The anelastic approximation relies on the hypothesis that the thermodynamics variables will not depart very far from a "reference state" defined as an atmosphere at rest in hydrostatic equilibrium. A variable V is therefore decomposed as V (x, y, z, t) = V ref (z) + V ′ (x, y, z, t), where V ref is the reference state's value of the variable and V ′ , the perturbation from this state. Meso-NH uses the Exner function Π, defined as Π = (P/P ref ) (R/Cp ) where P is the pressure, P ref = 10 5 Pa is the reference state's pressure, R the gaz constant for dry air and C p the specific heat content pressure for dry air. θ, the potential temperature, is 10 defined as θ = T /Π where T is the air temperature.
Under the hypothesis made in our study, the model's equations are written as follows: Where ρ (kg m −3 ) is the air density, θ (K) and Π (J kg −1 K −1 ), the potential temperature and the Exner function previously  (4) the thermodynamic equation.
The following two subsections describe the implementations performed in Meso-NH to make it a suitable tool to simulate flows generated in laboratory experiments -namely, the addition of the viscous terms in the model's dynamic, and the implementation of an explicit no-slip bottom boundary condition.

Addition of the viscous diffusion terms
According to the Kolmogorov theory of turbulence (see Kolmogorov, 1941a, b), the turbulence kinetic energy is transfered from the largest turbulent eddies to the smallest ones until it is dissipated by viscosity. This transfer is called the energy cascade. Kolmogorov (1941a, b) postulated that when the flow is turbulent enough, the scale of viscous dissipation -called 30 the Kolmogorov microscale, l d (m) -depended only on the the kinematic viscosity of the fluid ν (m 2 s −1 ) and ǫ, the spectral rate of turbulence kinetic energy transfer, which is equal to the production rate of turbulence kinetic energy. The Kolmogorov microscale l d can be expressed as: where l e is the scale of the largest turbulent eddies and Re is the dimensionless Reynolds number. Re is defined as where U (m s −1 ) is the flow's velocity, Λ (m) is a typical length of the flow and ν (m 2 s −1 ), the kinematic viscosity. Re rep-5 resents the ratio of the intertial forces to viscous forces, it is commonly used to distinguish laminar and turbulent regimeshigh Reynolds numbers being associated with turbulent regimes. In typical turbulent atmospheric flows, the Reynolds number reaches values of 10 7 and above. The associated Kolmogorov microscale l d is then equal to 1 mm or less.
In meteorological models, the viscous dissipation is not taken into account. Given their typical resolution, they can not represent this phenomenon. At best (in LES), part of the turbulence remains subgrid and has to be parametrized. The dissipation of energy at the end of the spectrum is then ensured by a numerical diffusion scheme. But in our case, Meso-NH is meant to be applied to laboratory flows whose Reynolds numbers are smaller than atmospheric ones, mostly because of the reduced dimensions of laboratory installations. Typical values of Re in experimental flows range from 10 2 to 10 4 (e.g. Eiff and Bonneton, 2000;Knigge et al., 2010). Baines and Manins (1989) have suggested that Reynolds numbers exceeding several hundred are necessary for the flow generated in laboratory to be similar to the atmospheric one. So even though Reynolds numbers in lab-15 oratory experiments may be lower than those observed in the atmosphere, the two flows can be in a similar regime (Townsend, 1980;Wyngaard, 2010) -this is in fact closely related to the processes one studies and the critical value may depend on other parameters (see for example Brethouwer et al., 2007). In such cases, the experimental Re values, associated with the smaller dimensions of laboratory, make it conceivable to explicitly resolve the turbulence up until the Kolmogorov microscale with a numerical model, which would also represent the viscous diffusive terms of dissipation.

20
To do so in Meso-NH, we added the viscous diffusion terms of momemtum F m and heat F t to the momentum equation (3) and the heat equation (4), following the Navier-Stokes equations: where P r is the Prandtl number, defined as the ratio of the momentum diffusivity to the thermal diffusivity.

25
As explained in introduction, the reproduction of atmospheric flows in laboratories rests upon the concept of dimensional analysis. Conversely, an experimental flow generated in a water flume or a wind tunnel can be transposed into an equivalent atmospheric flow, provided that the dimensionless numbers characterising the flow -such as Re -remain the same in both cases. So when simulating a laboratory experiment in Meso-NH, we want to have exactly the same Reynolds number in the 30 numerical simulation. This point raises no difficulty when the model is run directly at the scale of the laboratory, but it does when working with the dimensions of the atmosphere where the typical Reynolds number are higher. To obtain a reduced Reynolds number, one can increase the viscosity ν in the model and keep typical atmospheric values for U and Λ. In the simulated atmospheric flow, the reduced Re will translate in to a bigger Kolmogorov microscale, making it also conceivable to run a DNS at the dimensions of the atmosphere.
To give an example, let us consider an atmospheric flow characterised by a mean velocity of U = 10 m s −1 and a typical length of Λ = 1000 m corresponding to the size l e of the largest turbulent eddies. The viscosity of the air mainly depends on the temperature ; for T = 15°C, ν = 1.5 10 −5 m 2 s −1 . The Reynolds number is approximately equal to 10 9 for this flow, and the a typical length of 1000 m, the viscosity has to be equal to 1 m 2 s −1 , and the Kolmogorov microscale is now equal to 1 m, making it also possible to consider a DNS.

15
In laboratory experiments simulating atmospheric flows in wind tunnel or water flumes, a boundary layer always develops above the ground and it interacts with the flow, as it does in the atmosphere. For example, it can have a significant effect on the orographic waves occurring in stably stratified flows. This interaction is the subject of numerous studies (e.g. Eiff and Bonneton, 2000;Vosper, 2004). The vertical development of the boundary layer can be enhanced with the use of elements of rugosity in the experiments and reduced with the use of a smooth surface. To isolate the effects of the boundary 20 layer, it would be necessary to suppress it in the experiments. But even with a very smooth surface, this is not possible. It can however be done with a numerical model, and this is one of the foreseen application intended to be performed with Meso-NH.
Anyhow, Meso-NH must be able to simulate the boundary layers forming in laboratory experiments.
The native bottom boundary condition in Meso-NH is a "free-slip" one: the normal component of the wind is set to zero at the ground level while no condition is applied to the tangential component of the wind. It corresponds to the following equation: Despite this free-slip condition, Meso-NH takes into account the ground friction through the surface wind shear u * , expressed as a function of the roughness length z 0 characterising the rugosity of the surface. This formulation requires to activate a turbulence scheme and to prescribe surface turbulent fluxes. In most -if not all -simulations run with Meso-NH, these basic requirements have raised no question so far. But in our case, Meso-NH is intended to be used explicitly to represent as 30 realistically as possible the boundary layer close to the bottom boundary condition. Therefore, the ground friction needs to be implemented explicitly with a no-slip bottom boundary condition, reading as follows: This bottom boundary condition would be trivial to implement in Meso-NH if the ground level was defined for the three components of the wind. But on the Arakawa-C grid (see figure 1) of the model only w, the vertical component of the wind, is defined for z = 0. u and v, the horizontal components of the wind, are located at the altitudes −∆z/2 and ∆z/2. Physically, the no-slip bottom condition impacts the rest of the flow through the advection of momentum and the viscous diffusion of the wind. In the free-slip condition, w is already set to zero for z = 0. To implement the rest of the no-slip condition, the velocity 5 of the exterior point (z = ∆z/2) is set equal and opposite of that of the interior value (z = −∆z/2): The no-slip condition can either be applied to the whole domain of integration or only to a subdomain. This feature can be necessary when for example, the flow is created by an obstacle towed in a laboratory experiment, the bottom friction is present only above the obstacle surface. Several ways of defining the subdomain on which the no-slip condition is to be applied have been implemented. The first one is based on a condition upon the height of the surface z s : the no-slip condition is applied 15 where z s > z 0 (see figure 1 for the definition of z s ). One can also define a subdomain with indexes in the x and y direction (rectangular subdomain), or with any given mask of points.
3 Validation against exact solutions of ideal flows 3.1 Validation of the viscous diffusion terms 20 We will only validate the viscous diffusion of momentum, the viscous thermal diffusion having been implemented in the same manner.

Ideal flow and analytical solution
The ideal flow we consider here is defined by the following initial state: where (x, y, z) refers to the cartesian coordinate system, θ 0 (K) and U 0 m s −1 ) are constants and L (m) is the length of the domain in the y-axis direction (the domain is infinite along the x-axis and the height z does not have to be specified for this Neglecting the Coriolis force, the momentum equation is simplified into: With a free-slip bottom boundary condition, this equation has the following exact solution:

Model set-up
Meso-NH is run with dry air, with the Boussinesq approximation and without the Coriolis force and no parametrization. The initial velocity of the wind, U 0 , is taken equal to 10 m s −1 and the kinematic viscosity is set at 100 m 2 s −1 . With this high viscosity, the extinction term e −νπ 2 L 2 t in the expression of u reaches a significant value within a fairly short period of integration. 10 In the y-direction, the lenght L of the domain is 4800 m, with rigid wall lateral boundary conditions. In the x-direction, the lateral boundary conditions are cyclic, which is equivalent to an infinite domain of integration along x. As the dimension in this direction does not matter, it consists of only 8 grid points. Since the flow is steady along the z-axis, 8 vertical levels are defined. The upper boundary condition is modelled with a rigid roof and a free-slip condition. The resolution is equal to 50 meters in every direction, and the time step is set to 5 s.

Results
Simulated values of u (at the first vertical level) are compared with analytical solution on figure 2. u is represented as a function of y for t = 0, and after 1, 2 and 3 hours of simulation. Results show that the model is very close to the analytical solution.
For the first 2 hours of simulations, the curves superimpose so well that one can barely see the difference. After 3 hours of simulation, the difference is more visible, but the relative error stays below 2.7%.

20
From this validation, we conclude that our implementation of the viscous terms is correct.

Ideal flow and analytical solution
This flow is known as the "Stokes's first problem" (Schlichting, 1968). It occurs over a suddenly accelerated flat plate. Over time, a boundary layer develops above the plate and it can be analytically described.

25
With an infinite bottom surface started impulsively from rest to a speed of U 0 in the x-axis direction, we obtain the simplified momentum equation: with the following boundary conditions: If we define its height δ as the height for which u(δ) = 0.01, we obtain: See Schlichting (1968) for a demonstration of this result.
This problem is equivalent to the one where the fluid would be set at the uniform speed of U 0 , with a free-slip on the bottom surface, and then suddenly submitted to the ground friction. In this case, the momentum equation is the same as (16) and the boundary conditions become: With a boundary layer thickness δ defined as u(δ) = 0.99 U 0 , we still have δ ≈ 3.64 √ νt.
In Meso-NH, the ground can not move. But this sudden shift from a free-slip to a no-slip bottom condition is exactly what happens when the model starts from an initial state where the velocity is set at a uniform value over the whole domain of integration.

Model set-up
The initial wind velocity, U 0 , is taken equal to 10 m s −1 . Three different values of viscosity are considered: These values of viscosities match the range of values obtained for typical Reynolds numbers in laboratory experiments (ν 3 is equal to the kinematic viscosity that will be defined in section 4 for the validation against experimental data).
Here again, Meso-NH is used explicitly. It is run for 1000 s, with the fourth-order centred scheme to transport momentum, inducing a time step of 0.5 s. Since the flat ground is infinite, the flow is essentially one dimensional. In the model, this feature is obtained with cyclic boundary conditions in both horizontal directions. Thus, the horizontal size of the domain does not 25 matter, and the resolution is not crucial either. 10 points are defined in each direction and the horizontal resolution is equal to 1 m. In the ideal flow, the velocity in the free stream above the boundary layer remains equal to U 0 . In the model however, the decrease of velocity in the boundary layer will be compensated by an increase in the upper levels for conservation purposes.
To limit the impact of this effect on the comparison to the ideal flow, the domain of integration must be high enough. It is set to 2000 m. The vertical resolution is equal to 1 m in the first hundred meters where, according to equation (19), the boundary layer is expected to develop during a 1000-second long simulation. Above, the resolution ranges from 2 m to 10 m, and it is equal to 50 m in the upper levels.

Results
The results are shown on figure 3, where the simulated and analytical values of δ, the boundary layer height, are plotted against 5 t, the time of integration. In all three simulations, the model's values of δ stay very close to the ones analytically computed.
The differences are slighty higher in the third simulation, where the kinematic viscosity is the smallest (ν = 0.064 m 2 s −1 ), which could indicate that the resolution of 1 m is not quite high enough for this case. But even so, the solution provided by the model stays very close to the analytical one, with relative differences smaller than 5 %. This facility is unique in particular due to its capacity to generate high-Reynolds numbers stratified flows with low confinement 15 effects. The range of heights of the boundary layers typically developping in this facility allow the use of reduced size models with scales ranging from 1 : 50 to 1 : 10000.
The laboratory also includes two smaller water tanks (4-meter and 7-meter long) as well as a 2.5 meters rotating turntable. The research activities of the laboratory ended at the end of 2014 due to severe budget cut, but part of the equipment is kept in use for lectures. has been chosen in order to isolate the effects of the ground friction from other processes and make a proper evaluation of the numerical code. In this rather old experiment (Perrier and Butet, 1988), the CNRM large stratified water flume was used with a 1-meter deep single layer of water (no density stratification) circulating at a velocity of 25 cm s −1 . The relative mean velocity 25 fluctuation was less than 1 %, thanks to the controlled stability of the pump regime (1/1000) and for the water temperature (less than 0.1°C per day).
The first five meters of the flume (on the right on figure 6) are used to laminarize the flow coming out of the pumps and pipes located under the canal. It is composed of three grids with a decreasing mesh size going from 20 mm to 0.5 mm. At the outflow of this 5-meter chamber, the turbulence intensity is less than 2 % -to be compared to the 20 % measured at the entrance of the The vertical profiles of u show an increase of velocities at the vicinity of the floor, with a maximum u max superior to u ∞ , the mean velocity above the boundary layer (see figure 7). This is due to the experimental configuration where the floor is actually 10 a false floor. This feature of u(z) is more pronounced for the profiles close to the entrance of the flow and it almost disappears for x ≥ 9 m. It makes it difficult to define a boundary layer height that would be consistent along the x-axis. In the previous sections, we used a boundary layer thickness δ 0 defined as: But in the present laboratory experiment, when the vertical profiles of u show a significant increase of velocity near the 15 floor, this classic definition does not make sens anymore. To overcome this problem, two other criteria were defined by Perrier and Butet (1988). These are δ 1 and δ 2 , defined as: For x ≥ 7.5 m, we were able to define δ 0 as u(δ 0 ) = 0.99u ∞ by considering the values of u for z > δ 2 . These three definitions 20 of δ are illustrated on figure 7. Wherever δ 0 can not be defined because of the false floor effect, it seems reasonable to assume that if this effect was suppressed, the value of δ 0 would lie somewhere between δ 1 and δ 2 . Therefore, we will use δ 0 = 0.5(δ 1 + δ 2 ) in these cases.

Model set-up and numerical simulations
As explained in introduction, a flow generated in a hydraulic flume can be transposed into an analog atmospheric flow. To -R e = U.Λ/ν As long as these numbers have the same value in the atmospheric flow and the experimental one, both problems are equivalent. the laboratory, all results will be presented accordingly. We just provide the corresponding atmospheric dimensions because they are more meaningful to atmosphere modellers (see table 2).
Our goal is to simulate the experimental flow with a DNS (Direct Numerical Simulation) that would accurately represent the observed turbulence, without any parametrization. We performed such a simulation -that we will simply call DNS -with a horizontal resolution of 1 mm and a vertical resolution ranging from 1 mm in the boundary layer to 10 mm at the top. The  it is not possible to run Meso-NH at this resolution. However according to Yakovenko et al. (2011), it is possible to run DNS for resolutions which are around 10 times the Kolmogorov microscale. We argue -and results will confirm -that the resolution of 1 mm, which is 20 times the Kolmogorov microscale, is sufficient.

5
An overview of the simulated flow is presented on figures 8 and 9 with a vertical slice of the velocitie U along the (x) direction, and several three-dimensional views of U . They were computed with instantaneous outputs at the end of the simulations (t = 82 s). These plots show the thickening of the boundary layer along x and the turbulent nature of the simulated flow. They reveal small-scale 3D isotropic patterns of turbulence in the boundary layer. This indicates that our version of Meso-NH might indeed be able to resolve the whole spectrum of turbulence for this flow. To determine if this is actually the case, we will now 10 compare the DNS results to the laboratory measurements in terms of boundary layer heights, vertical profiles of velocity and turbulence intensities, and other characteristics of the turbulence.
To compute further diagnostics in DNS, instantaneous velocities will need to be averaged over time in stationary regime.
A statistical test is performed to determine the time t s when the flow becomes stationary. We consider u the x-component of this parameter, we also plotted the observed values of δ 1 and δ 2 , as well as the simulated δ 1 (there is no δ 2 in the simulation).
DNS stays fairly close to the experimental data. It follows almost exactly the observed values in the first 2 meters. Up until x = 6 m, the simulated boundary layer stays within the range of uncertainty given by δ 1 and δ 2 , after which the model tends to underestimate the height of the boundary layer. Overall, with errors staying below 18 % on δ 0 , one can say that DNS provides rather realistic values of the boundary layer height.

30
We will now examine vertical profiles of u and turbulence intensities of DNS for a further assessment. Figure  Along with mean velocity profiles, the turbulence intensity I (I = √ u ′2 /u) was also measured in the laboratory experiments.
In the DNS simulation, these quantities were directly computed with the instantaneous outputs and the values of velocity averaged over a 5-second time period. The comparison is presentend on figure 13. It shows that the simulated profiles of I are 10 quite realistic in the boundary layer. Above, the turbulence of DNS is virtually null. This is not surprising, as the two sources of turbulence are the white noise applied on the inflow (for x = 0 m) and the vertical sheer of velocity which is present only in the boundary layer.
The von Karman and Prandtl theory gives a universal law for u in the surface boundary layer in neutral conditions of 15 stratification: with κ, the universal constant of Karman κ ≈ 0.41 ; u * = √ u ′ w ′ , the mean turbulent flux and z 0 , the roughness length characterising the rugosity of the ground. In situ observations in the atmosphere follow this universal profile. It is valid for z >> z 0 and z << δ 0 , that is to say not too close to the ground or the top of the boundary layer.

25
In the model where the ground is perfectly smooth, it is not possible to give a theoretical value of z 0 . However, a logarithmic regression can be performed on the vertical profiles of velocity, from which estimates of u * and z 0 can be deduced -they will be refered to as u * and z 0 . We did so in DNS and for the experimental data, for values of z comprised between δ 0 /10 and is disturbed by the false floor effect but for these two profiles, the logarithmic regression does not make much sense anyway.
Elsewhere, u * tends to be slightly higher than u * , in both DNS and the data. This means that the simulated and experimental profiles of velocity in the boundary layer are steeper than what the von Karman and Prandtl theory predicts. Regarding z 0 , the values deduced from the logarithmic regression have the same order of magnitude in the experimental data and in DNS (again, the first two profiles are discarded because of the false floor effect).

5
To sum up, results showed that Meso-NH is able to accurately simulate the turbulent flow observed in the laboratory. With the implementations presented in this study, Meso-NH has thus proven to be successful in explicitly representing the turbulence developing in laboratory experiments. With a computational domain of half a billion grid points, this DNS was made possible thanks to the adaptation of Meso-NH for massively parallel computing (Pantillon et al., 2011). Nonetheless, the DNS simulation we performed in this study can open some new prospects in the development of Meso-NH for the future. For instance, equivalent DNS simulations could be used as references to assess the parametrization of the turbulence based on the prognostic equation for the turbulence kinetic energy (TKE) with a closure by mixing length (Cuxart et al., 2000) in Meso-NH. A new mixing length (Rodier et al., 2017) has been recently implemented in the model for neutral stratification. It could be evaluated with the DNS in terms of turbulent intensity and eddies size. Hence, a DNS configuration could prove very 5 useful to improve and validate turbulence schemes in LES, the same way LES simulations have been used for years to develop and improve convection paramatrizations in GCM (General Circulation Models) and NWP (Numerical Weather Prediction) (see for example Siebesma et al., 2003;Rio and Hourdin, 2008). In any case, we believe the possibility to run Meso-NH in a DNS mode constitutes a noteworthy advance in the model's development.
Code availability. The model is freely available and can be downloaded at http://mesonh.aero.obs-mip.fr/mesonh53/Download       Table 4. Values of u * and z 0 deduced from the logarithmic regressions of the observed and simulated profiles of mean velocity u(z) (figure 14).