SURFACE BOUNDARY LAYERS THROUGH A SCALAR EQUATION WITH AN EDDY VISCOSITY VANISHING AT THE GROUND

. We introduce a scalar elliptic equation defined on a boundary layer given by Π 2 × [0 , 𝑧 top ], where Π 2 is a two dimensional torus, with an eddy vertical viscosity of order 𝑧 𝛼 , 𝛼 ∈ [0 , 1], an homogeneous boundary condition at 𝑧 = 0, and a Robin condition at 𝑧 = 𝑧 top . We show the existence of weak solutions to this boundary problem, distinguishing the cases 0 ≤ 𝛼 < 1 and 𝛼 = 1. Then we carry out several numerical simulations, showing the ability of our model to accurately reproduce profiles close to those predicted by the Monin–Obukhov theory, by calculating stabilizing functions.


Introduction
This paper is devoted to study a scalar elliptic equation which parameterize the mean velocity of the air in the atmospheric Surface Boundary Layer (SBL), where the turbulent stresses are balanced with the friction forces on the ground.This is part of the more general framework of the turbulent boundary layers, initially developed by Prandtl [19], then by von Kármán [26], who highlighted the role of logarithmic profiles relative to the height in such layers (see also in [13,24,25]), called the log-law, which was validated by several numerical simulations, for instance by using turbulence models such as the  −  model (see [18] and further references inside) and-or by stochastic models [17].
The Monin-Obukhov theory [14] states that under non-neutral conditions, the mean velocity profile differs slightly from the log-law, the difference being determined by stabilization functions.This theory is used in much more general (SBL) regimes [16], and is the basis of most atmospheric flow simulations near the ground, which raises the question of the determination of the stabilization functions.
The starting point is the 1D differential equation that yields the log-law from a theoretical point of view [24], namely together with appropriate boundary conditions, where  = () denotes the mean horizontal velocity component that is supposed to only depend on the height  > 0; in this framework,  ⋆ is the friction velocity and  the von Kármán constant.We wonder if a similar simple multi dimensional PDE satisfied by  = (x ℎ , ), x ℎ ∈ IR  ( = 1, 2), is able to yield Monin-Obukhov profiles type, which suggests to introduce the equation, where  is the Boussinesq force specified by a temperature supposed to be constant in this paper,  turb =  turb () is an eddy viscosity and  ℎ > 0 an horizontal viscosity coefficient.According to standard assumptions and dimensional analysis [5,10,15], we should have in the domain 0 <  <  top , where  top denotes the height of the (SBL).However, we know that eddy viscosities that vanish at the boundary are source of serious mathematical issues [2][3][4] and are often studied by means of weighted Sobolev spaces [8].Moreover, the case given by (1.2) is critical as we will see later in this article, giving very weak solutions with only  1/2 regularity, which does not allow to set a boundary condition at  = 0.This is why in many models,  turb is taken to be constant in a viscous sublayer [0,  0 ].The same issue occurs in the case of the Smagorinsky model [23], where the eddy viscosity denoted by  smag is given by  smag =  2 |  | near the ground  = 0.This is why in the Smagorinsky case, several authors have suggested to replace the physical  smag by  smag =  2(1−) 0  2 |  | for some 0 <  < 1 [4,[20][21][22], to obtain more regular solutions and to be able to take into account appropriate boundary conditions.This suggests to consider in our case general eddy viscosities of the form There is the question of the boundary conditions.It is natural to set  = 0 at the ground  = 0. Following [9], we take a Navier friction condition at the top of the boundary layer  =  top (also named Robin law), which is a fairly transparent condition, easy to deal with in a variational formulation.In order to be consistent with the numerical simulations, we take periodic boundary conditions in the horizontal directions.Therefore, the PDE problem we consider in this paper is the following 1 , where x ℎ = (, ) and 0 <  ≤ 1, for a given 2D torus denoted by Π 2 , where the term ,  ≥ 0, is a stabilizing term, useful only in the case  = 1,  > 0 and   > 0 are given coefficients that will be calibrated by numerical experiments.
In this paper we prove the existence of a weak solution to problem (1.4) in an appropriate weighted space for 0 ≤  < 1, see Theorem 3.1 below.Then for  = 1, we prove the existence of a very weak solution  ∈  1/2 (BL) for  > 0, in Theorem 4.1, which is the main result of this paper.In this result, we do not impose  = 0 at  = 0.It is based on a Neças lemma type (also known as Lions Lemma), Lemma 4.1.In this case, the difficulty is due to the mixed boundary condition and we cannot directly apply the results of [1,6].We had to make many efforts to prove this essential result in this problem, based on the interpolation theorem [12].Finally we carry out several numerical simulations based on a Freefem++ code [7], which allows to evaluate the difference between the solution and the log-law.In particular we observe that, despite the lack of theoretical regularity, the physical case  = 1 remains the most accurate to parameterize the SBL, and we are able to calculate numerically in several different regimes the stabilizing functions, given by formula (5.8) below, which validates model (1.4) in terms of the Monin-Obukhov theory.
The paper is organized as follows.In a first part we develop the modeling sketched in the introduction, and we set the physical constants involved in the simulation.Then we study the case 0 ≤  < 1, by viscous regularization and proving that the natural weighted space related to the problem is embedded in a standard Sobolev space  1, (BL).Then we focus on the case  = 1.A large part of the paper is devoted to the study of the function spaces and Neças Lemma by means of Fourier series, which allows to prove that the natural weighted space related to the problem is embedded in  1/2 (BL).The last section of the paper is devoted to the numerical results.

Modeling
This section is devoted to recall some basic elements of the theory of turbulent boundary layers, and to fix the general framework of the model which one studies.The steady mean fluid velocity in such boundary layer, denoted (BL), is denoted by instead of u for the simplicity, where u ℎ = (, ),  ∈ ]0,  top [,  top > 0 being the bottom of (BL).For instance, if (BL) models the surface boundary layer,  top ≈ 100 m.We also will need to consider the roughness length  0 , which depends on the nature of the ground, and varies from 0.0002 m in open sea, to 1 m for city centre with high-and low-rise buildings.Note that we are in a flat domain and the splitting of both variables and unknowns into horizontal and vertical will be of particular use to identify the problem and give a clear formulation.

Assumptions, general equation and issues
Let  > 0 denotes the kinematic viscosity of the fluid.It is commonly supposed that in standard (statistical) steady BL it hold the following: -the pressure is constant and the vertical part of the mean velocity vanishes, that is  = 0 and, even if it means making a change of coordinates, we can assume  = 0; -the mean velocity u = (, 0, 0) depends only on the altitude, that is  = (); -the eddy viscosity  turb depends on  and  ⋆ = √︀ |  (0)|, the so-called frictional velocity, which is the tuning parameter of the system (see [18,24]), which yields where   ≈ 15,  ⋆ ≈ 10 are non dimensional constants, that we have calibrated by numerical simulations.Typical values of  ⋆ range from 2 to 10 ms −1 .-all terms in the fluid equation are negligible compared to the turbulent diffusion term.
These assumptions lead to the following equation for the mean velocity u = (, 0, 0), which formula, once integrated between a given  0 and  top with appropriate boundary conditions, yields the well known log-law, uniform in x ℎ , using the calibration constants   and  ⋆ : Generally, for  ∈ [0,  0 ], called the viscous sub-layer, a linear profile is considered such that  = () is continuous over ]0,  top [, and (0) = 0, which means Let  Log denotes the function defined over ]0,  top [ by (2.3) and (2.4).When the stability of the atmosphere is non-neutral and due to the effect of convection, which means that there is a non zero source term in (2.2), stabilizing functions must added to  Log to get the right velocity profile, according to the Monin-Obukhov theory [14], which means that where the function Ψ(x ℎ , ) is deduced from similarity arguments or from experimental data.Examples of such stabilizing functions can be found in [16].
Remark 2.1.Normally in usual industrial models,   stands for the von Kármán constant, the value of which being equal to 0.4, and  ⋆ = 1.However, due to the scales of our simulations, we have to take other values of these constants to get numerical results related to the physical data in the atmospheric SBL.
Our aim is to find a comprehensive PDE model, such that: (1) is defined over ]0,  top [; (2) includes an eddy viscosity of the same form as that given by (2.1), where the profile  = (x ℎ , 0) also depends on the horizontal variable; (3) is able to calculate stabilizing functions such as in (2.5) in various atmospheric regimes.
Before embarking on nonlinear complicated 3D equations of fluid mechanics, we consider as a first step the following elliptic toy-model in BL = IR 2 × (0,  top ): for some  > 0,  ℎ > 0, ∆ ℎ =  2  2 +  2  2 .The term  in (2.6) stands for a numerical artefact of an evolutionary term   , and serves as a system stabilizer, especially in the case  = 1.It can be taken equal to zero in the finite element simulation thanks to the numerical dissipation, due to the discretization.
In physical applications, the source term  is the Boussinesq force, namely where  is the temperature of the fluid,  0 its value at the ground,  ≈ 10 ms −1 is the gravity coefficient,  is the coefficient of thermal expansion, a typical value of which for dry air is varies between 0.002 K −1 and 0.003 K −1 .

Boundary conditions
The choice of the boundary conditions may be an issue, and there are many options.We consider the case of a BL that flows over a rigid wall, which means that we take an homogeneous boundary condition on the bottom Γ  , (x ℎ , 0) = 0.Moreover, we consider that this BL is coupled on top with a layer of fluid which exerts a frictional force on it.Therefore, as in [9] one can can take a linear Navier-Boundary condition like where   > 0 is a frictional coefficient and  =  (x ℎ ) is the velocity of the top layer.In the numerical simulations we have taken , where  Log is that given by (2.3) and (2.4).The coefficient   will be numerically optimized in function of  ⋆ , (x ℎ ) is a small perturbation term.
Remark 2.2.According to the results of [9], we expect that for large values of   , the boundary condition (2.8) converges (in some sense) to the continuity condition (x ℎ ) =  (x ℎ ) at Γ top , which is well confirmed in this framework by the numerical simulations.

Alternatives and general framework
As we will see in the following, the eddy viscosity given by (2.1) yields variational (or weak) solutions to problem (2.6) that are in  1/2 , and not much more.In particular the homogeneous Dirichlet boundary condition at the bottom cannot be checked, which is consistent with the log-law.This is why we ask the question whether or not it is possible to identify alternate eddy viscosities, close to (2.1) but giving more regularity to the system, being not critical for the notion of trace, and so on.We wonder if that gives good approximations of the usual BL profiles.In this way, it is natural to consider  turb of the form (2.9) the main feature of which is that it degenerates at the ground, but with a different velocity.The parameter  0 has the dimensions of a length and it is needed to have a consistent expression for the viscosity.We take as boundary conditions on the bottom and on the top, we write the friction law (2.8)like a standard Robin condition in the form: for  =  1− 0    ⋆ and where  =    .It remains to clarify the boundary conditions in the x ℎ -axis.For practical numerical simulations, we have to limit ourselves to a finite computational box [0,   ] × [0,   ]×]0,  top [ for given scales   and   , which raises the question of the boundary conditions at the entrance, exit and sides of the computational box, namely In [11] we have considered at Γ in a fixed given field in  1/2 00 (Γ in ) and nonlinear Neumann transparent boundary conditions at Γ 0 which yields serious technical issues, both theoretically and numerically.In this paper, we opt for horizontal periodic boundary conditions, which means that  must satisfy for all formal derivative operator   ,  ≥ 0,   ( +   ,  +   , ) =   (, , ), and imposing the invariance for translations implies working with the torus Π 2 given by where  2 denotes the set of wave vectors given by: This setting is usual in practical numerical simulations for technical convenience but not only since also it eliminates the problem of an infinite domain in the analytical studies.For instance, it was used in [18] for simulating a boundary layer by a RANS turbulent model.

The 0 < 𝛼 < 1 case
Throughout this section, we assume that 0 <  < 1 and  > 0 is fixed and we study the problem in a simpler setting.In fact, in this case one can reduce the problem to another one in a setting of unweighted Sobolev spaces for which both existence and interpretation of the solution are standard.Later on we will see the possible treatment of the limiting case.
We will use standard Lebesgue and Sobolev spaces and from now on We specify the function spaces we are working with, and then we prove the existence of a solution by a viscous regularization.To start with, we observe that any strong solution to problem (1.4) satisfies the energy balance:

Function spaces
In this section, we define the function spaces to be utilized throughout the paper.We employ the conventional notations within the framework of distributions, adapted to accommodate the mixed boundary conditions under consideration.
-Let   (BL) denote the space of functions  ∈  ∞ (BL) such that there exists . This space is endowed with the standard Schwartz topology, characterized by the semi-norms ‖  ‖ ∞ .In simpler terms,   (BL) consists of  ∞ functions that vanish in a neighborhood of the bottom Γ  of the boundary layer.-Let  ′  (BL) denote the topological dual of   (BL), a sort of distributional space that inspires our chosen notation.
- 1, 0, (BL) ( > 1) denotes the closure of   (BL) for the norm ‖∇ ℎ ‖ 0, + ‖  ‖ 0, , where ‖ • ‖ , stands for the usual  , norm.In particular, According to the energy balance (3.2), we are led to consider the natural weighted space   defined as the closure of   (BL) equipped with the norm and the following inequality holds ∀  ∈   , for some constant  > 0.
Proof.Let  > 0 and  > 0 that will be fixed later.Let  ∈  0, (BL).The Hölder inequality yields: The second integral is well-defined if and only if 2 2− < 1.Then, choosing  such that 2  =  yields the condition  < 2 +1 .Inequality (3.5) follows after an elementary calculation and integration with respect to the dx ℎ variables.
It follows from Proposition 3.1 and standard reasoning on Sobolev spaces, that functions in   have a trace at  = 0 equal to zero, and also we have the following characterization

Weak formulation
Proposition 3.1 can be rephrased, to work in standard (unweighted) Sobolev spaces, as follows which is put in duality with the set Throughout the section, we assume that ) The following definition of a weak solution to problem (1.4) is motivated by standard rules about integration by parts, combined with the boundary conditions under consideration.
Remark 3.1.Note that all the terms in the integrals written in (3.11) are well-defined.However, the solution  cannot be a priori taken as test function, which is an issue.As a consequence, we are not able to prove the uniqueness of this solution, even if the problem is linear.
Before all, we notice that combining the energy balance (3.2) and (3.5) with standard calculus inequalities, yields for any 1 <  < 2 1+ the following estimate in  1, 0, (BL), satisfied by any given regular solution  to the variational problem (3.11): where   → ∞ as  → 2 1+ .The aim of the rest of this section is proving the following existence result.Theorem 3.1.Problem (1.4) admits a weak solution  ∈   .Moreover, the solution satisfies the energy inequality

Viscous regularization and proof of the existence result
From now we set since this value is the critical one for the embedding of weighted Sobolev spaces.As a technical tool to prove existence of weak solutions, we regularize problem (1.4) by adding a viscous term in the -direction, which means that we consider the following problem, for a given  > 0: The existence and uniqueness of a weak solution to problem (3.16) is straightforward by the Lax-Milgram theorem.Moreover, as   can be taken as test function, it satisfies the following energy balance (equality): From this, we are able to finish the proof of Theorem 3.1 by taking the limit in (3.16) when  → 0. We deduce from (3.17) and standard calculus inequalities that the family (  ) >0 is uniformly bounded in   , as well as in  1, 0, (BL) for any 1 <  <  ⋆ say Arguing as in Chapter 7 of [24], we can extract a (sub)sequence (  ) ∈I N that weakly converges to some  ∈   , which is also weakly converging in  1, 0, (BL) for all 1 <  <  ⋆ , and which is strongly converging in  2 (BL).Moreover, by the trace theorem and the Sobolev theorem, the sequence (  ) can be chosen such that in addition ([  ]) ∈I N strongly converges to [] in  2 (Γ top ).Let 1 <  <  ⋆ and take as test function in formulation (3.16).We have to take the limit in the various terms of (3.16), which we do step by step, starting with the diffusion term.
Let  ∈   , and let the linear form Ψ  given by By the Hölder inequality we obtain hence the energy balance (3.13), which concludes the proof.
Remark 3.5.The solution we exhibit is obtained by viscous regularization.We also can think to directly get a solution by applying the Lax-Milgram theorem in the space   , which requires a different approach that we have voluntarily skipped here.As already stressed in Remark 3.1, we do not know if they are completely equivalent; we conjecture that this is the case.

The case 𝛼 = 1
We now consider the system: As we shall see, the best we can get is an estimate in an  1/2 -like space near the bottom Γ  .Therefore, we are not able to define the trace of such functions at Γ  .This is why we do not impose any boundary condition there, contrary to what we did for the case 0 <  < 1.Nevertheless, we still impose  |Γ  = 0 in the numerical code, since the finite element formulation yields additional numerical dissipation that sufficiently regularizes the system.
In order to derive this  1/2 estimate, in the same spirit as in [2], we first need the following version of Neças lemma, the proof of which will be complete by the end of this section.To define the function spaces we will utilize, let (BL) denote the space of functions  ∈  ∞ (BL) such that there exists  = () ∈ ]0, 1[, ensuring that the support of  is contained within Π 2 × [(), 1 − ()].This space is endowed with the standard Schwartz topology, characterized by the semi-norms ‖  ‖ ∞ .In simpler terms, (BL) consists of  ∞ functions that vanish in the vicinity of both the bottom and the top of the boundary layer.We use  ′ (BL) to denote the corresponding space of distributions.Additionally, so far no risk of confusion arises, we still set and  −1 (BL) ⊂  ′ (BL) is the topological dual of  1 0 (BL).
Corollary 4.1.Let  1 denotes the space: Then The space  1 is the space that naturally suits problem (4.1) according to the energy balance (3.2) that still holds when  = 1.It must be stressed that: -showing that any solution of problem (4.1) is well-defined in  1 by using (3.2) requires  > 0; -we have for any  > 0, which makes consistent the boundary condition at Γ top .
Definition 4.1.We say that  ∈  1 is a very weak solution to (4.1) if for all  ∈  1 , The main result of this section is the following.
Theorem 4.1.Let  ∈  2 and  ∈  2 .Then, there exists a very weak solution  ∈  1 ⊂  1/2 (BL) of problem (4.1), which satisfies We argue by approximation by using Fourier series expansions which allows for some explicit and direct computations.This is the reason why, before proving Lemma 4.1 and Theorem 4.1, we need to prove a bunch of convergence results about Fourier series in this context of mixed boundary conditions, periodic in the x ℎ -axis but not in the -axis.This will be the purpose of the next sections.

Framework
Let  2 = 2  Z × 2  Z and let k = (  ,   ) ∈  2 any wave vector.In the following, we set Let  ∈  1 (BL) and let  k =  k (; ) denotes the "horizontal" Fourier's coefficient at the wave vector k, namely where  = √︁  2  +  2  .Let   be its partial sum of the Fourier series defined by

.11)
We will prove the following result.
We first notice that the following isometries hold: Therefore, from the standard Neças inequality (see [6]), we easily get the following result.Notice that the explicit computations allow us to prove that the right hand side bounds the left hand side with multiplicative coefficient equal to one, independently of  ∈ IN.
Lemma 4.3.Let   ∈  ′  (BL) such that   , ∇  ∈  −1  .Then   ∈  2 (BL) and one has It remains to pass to the limit in (4.13) when  → ∞ to prove Lemma 4.1.We first must prove Lemma 4.2, which will be done step by step in the following subsections.

𝐿 2 convergence
We prove in this section the following convergence result, which is well-known but we give a self-contained treatment of all the results in this section.Lemma 4.4.Let  ∈  2 (BL),   its Fourier's series expansion as given by (4.8).Then   →  in  2 (BL) as  → ∞.
Step 2: Convergence.We deduce from classical results that for almost all  ∈ ]0,  top [, Put another way: Moreover, by the previous step, giving which concludes the proof.We must pass to the limit in the r.h.s of (4.24).By (4.20), we have Therefore, as   →  weakly in  1 0 (BL), By using the same argument, we also have the following, where we note    =   .
Proof.Let  ∈  2 (BL) and   ∈  2 (BL) given by (4.22).Then, we know from Lemma 4.4 that for any  > 0 be fixed,     converges to   in  2 (BL).We write Then, since ‖  ‖ ≤ 1, we get The rest of the proof is straightforward.

𝐻 1 Convergence
We are now able to prove Lemma 4.2, stating the convergence of the sequence (  ) ∈I N to  in  1 (BL).We already know from Lemma 4.7 that Then, for a.e  ∈ ]0,  top [,   is given by

.27)
As the Lebesgue measure over BL is -finite, we deduce from the Lebesgue-Fubini theorem that for a.e  ∈ ]0,  top [, (•, ) ∈  1 (Π 2 ).Therefore, always a.e  ∈ ]0,  top [, and by standard results, Moreover, by the Young inequality, where ‖.‖ Π2;, denotes the usual norm of the Sobolev space  , (Π 2 ).Therefore, by the Lebesgue Theorem, we get and by an analysis similar to the previous one, we can conclude that in  2 (BL), which finishes the proof of Lemma 4.2.
We now need the definition of the space  Then in our case, equipped with the norm According to Lemma 4.8, we have the following.
Now following [12], let us consider the linear operator Λ = ∇ −1 :  2 →  −1 , defined such that with natural notations to denote the various scalar products.Then, the space [ −1 ,  2 ] 1/2 is the domain of Λ 1/2 , and one has, according to the above, In order to complete the interpolation tool box, let us consider  −1 as the closure of  ∞ (BL) subjected to the norm 0 is the closure with respect to the norm Note that  0 ≈  1 (BL).The question is the characterization of the interpolation space [ −1 ,  0 ] 1/2 .From (4.33), we have which allows to characterize the space [ −1 ,  0 ] 1/2 thanks to (4.34) by its norm given by:

Neças Lemma and consequences
We are now able to prove in this section Lemma 4.1, which we recall that it states the following inequality, satisfied by any distribution  ∈  ′ (BL) such that , ∇ ∈  −1 (BL), the dual of the space  1 0 (BL) as defined by (4.2).Then, we will prove Corollary 4.1, namely with continuous injection, where  1 is defined by (4.4), equipped with the norm .
In particular, being  1 a closed subspace of  1/2 (BL) turns out to be an Hilbert spaces, allowing to use all the machinery of complete vector spaces, which makes possible the proof of Theorem 4.1.
Proof of Lemma 4.1.Let  ∈  ′ (BL) such that , ∇ ∈  −1 , and   , ∇  ∈  −1  be given by where  −1  is defined by (4.11).We have, for all  ∈ IN, and note that (by explicit computation) it easily follows that the constants on the right-hand side are equal to 1, independently of  ∈ IN.Therefore, according to inequality (4.13), Then, the sequence (  ) ∈I N is bounded in  2 (BL), and we can extract a subsequence (still denoted by (  ) ∈I N ) that weakly converges to some ṽ, which is equal to  in  −1 (BL), therefore in  2 (BL), which yields  ∈  2 (BL) and Proof of (4.37).Note that a similar result was already obtained in [2], using a former result of [1] combined with an interpolation argument.Going back to basics, we give here a self-contained and combined proof, by using the Neças lemma and the interpolation theory.Indeed, From Neças inequality (4.36), we deduce that the injection  :  −1 →  2 is continuous.Moreover,  :  0 →  1 is also continuous.Therefore, by the interpolation theorem [12], also the restriction of the identity (denoted still by ) : is continuous.In particular, by (4.35), there exists  > 0, such that

.39)
Conclusion of the proof.Thanks to (4.39), in order to prove (4.37), all we have to do from now is proving that the following inclusion holds true with continuous injection.In the case relevant for our problem, we have Let  ∈  1 , and  ∈  1/2 00 (BL).Then, 00 (BL)) ′ , and In the same way, we have also Combining all the previous inequalities yields , which concludes the proof.
The functional setting for the proof of Theorem 4.1 is more sophisticated than that used in the non-limiting cases 0 ≤  < 1.Despite the proof being a rather standard application of the Lax-Milgram lemma, the choice of the underlying function spaces is obliged by the nature of the problem and the fact that  = 1 implies that the function spaces do not embed in any standard Sobolev space with trace at the boundary.For this point cf.Proposition 3.1, which we recall is false for  = 1 and a counter-examples is easily built by means of a double logarithmic function.This is why we resort to the space  1 and, despite the abstract simplicity of the result, the interpretation of the notion of the solution is of particular difficulty, since the solution satisfies a weak formulation which is not the same as the strong one.As we will discuss later, the obtained solution has problem in the interpretation of the value at  = 0, for which the functional setting is not proper.Nevertheless, we cannot change it, since it is determined by the equations themselves, so we have to extract the maximum of information possible from the solution.
Proof of Theorem 4.1.We are now ready for the proof of the main result of the paper, that is the proof of existence of weak solutions in the limiting case  = 1.In this case once we have the adapted functional setting we can apply Lax-Milgram in the space  1 with the weak formulation defined by: ∀ ∈  Moreover, as is standard for these problems the function  ∈  1 can be taken as test function, proving the energy equality Moreover, this can be used also to shows uniqueness of the weak solution.The drawback is the impossibility of controlling the trace.
The only missing point is to observe that we proved negative norm lemma for functions vanishing at { = 0} ∪ { =  top }, while now we have a friction law at the upper boundary.This can be easily overcome by using the fact that the conditions are of Neumann (Navier) type at the top.Hence, considering the space instead of  1 0 (BL) as in (4.2), we can easily convert to the case of a Dirichlet problem in a doubled domain Π × (0, 2 top ), with a reflection along the line  =  top .Then, the proof remains the same as in the previous case.Note that since we have  1 functions for which the trace is not well-defined, we cannot expect that our weak solution satisfies  = 0 at the bottom boundary.So in this case the weak and strong formulation are not giving the same result and the same have been noted, in a slightly different setting, by Rappaz and Rochat [21], for the von Kármán problem.They also noted as the trace evaluated numerically is strongly depending on the mesh-size, as is expected, and that the value at the boundary is not under control.

Numerical experiments
In this part, we aim to check if the model gives a good approximation of the Monin-Obukhov law (2.3),depending on the values of .We solve the problem (1.4)  where and (5. 3) The aim of this numerical study is to check if the numerical solution  of the problem (1.4) approaches the known log-law, in a sense we will develop below.
First we will explain in Section 5.3 how to get vertical velocities from the Freefem++ resolution of the problem (1.4) and the tools to compare it with the log-law  Log .
Then, we will discuss the influence of the different parameters ( in Sect.5.2),   and  ⋆ in the Section 5.3 and calculate the difference between the numerical solution and  Log , which allow us in Section 5.4 to derive an analytical formula of the stabilization function Ψ deduced from the numerical results by interpolation.
Finally, we will consider the influence on the -axis of a small perturbations , as involved in the wind at the top given by (5.2).Remark 5.1.We use in the code the command  = 0 at  = 0, and we take  = 0. Contrary to what the analysis predicts, the case  = 1 works very well, even at high resolutions.As Freefem++ is a finite element software, we think that the numerical simulation involves a numerical dissipation which sufficiently regularizes the equation, even in the case  = 1.However, we did not have studied yet this numerical aspect of the problem.

Parameters
Different parameters will have influence on the simulation: some will be fixed, and some will be specifically studied.The size of the box will always be in the following [0, ] × [0,  top ], where  = 1000 m and  top = 100 m.
-After several simulations to get velocities that can be measured in situ, it looks like that the best values for the calibration constants  * and   were  * = 10 and   = 15, as already mentioned in Section 2.1.
Examples of log profiles have been plotted in Figure 1 with these values,  ⋆ = 10,  top = 100, and respectively  0 = 1 and  0 = 10.-The height of the viscous sub-layer  0 in the following simulations will be very small compared to  top , with a ratio smaller than 0.01.Besides, the viscous sub-layer height  0 will taken equal to 0.1, which respects a ratio  0  top = 10 −3 < 0.01.
-The parameter  is chosen equal to 0.
-The source  will be considered as constant:  = 5.According to formula (2.7), this means that  ≈ −23 ∘ C for  = 2 × 10 −3 , for instance in the case where the ground is at 0 ∘ C and the air is dry and and cold at a constant temperature equal to −23 ∘ C, a situation that can happen in the mountains.-The perturbation  which appears in the expression (5.2) of  will be taken equal to 0 in the next subsections, except in Section 5.5 (see (5.9)).
The velocity constant  * can be seen as a "wind regime switch" belonging to the speed range [2 ms −1 , 10 ms −1 ], which corresponds to what is generally measured for flows over rough grounds.This is the main parameter of our simulations, the influence of which will be studied in the Section 5.3.
-Finally, the frictional coefficient   will be chosen according to the  ⋆ .In the Section 5.3 we will show that   ≃ 10 6 is giving convincing results.

Errors
Let   denotes the vertical velocities at   = , given by, for every  ∈ [0,  top ]: which will be more relevant than the error because of the importance of the velocities, when  ⋆ > 7 ms −1 for instance.
Without the perturbation , the vertical velocities   are very close each other, as we can see in Figure 2 and in the Table 1.Even if the difference between the vertical profiles is small, we consider the mean vertical velocity and we will use it instead of the   .To perform the simulations,  = 11 and   will belongs to {0, 100, 200, . . ., 1000}.

Influence of alpha
We observe that the more  is getting close to 1, the better the model is, as we can see in Figure 3: this is correct in the sense that the difference between the calculated profile and the log profile is smaller.The blue curves correspond here to the values (), where  ∈ {1, 2, . . ., 100}.
As a result, the model is relevant only for  = 1.We will take  = 1 in the next simulations.

The three different regimes: influence of 𝐶 𝐷
We have calibrated the constant values to get wind velocities which are physically relevant (in ms −1 ).We have observed three different regimes for  * respectively equal to 4 ms −1 , 7 ms −1 and 10 ms −1 corresponding to small wind, medium wind, and storm wind.We consider for each case the errors and relative errors we get in function of the   coefficients.We can see in Figure 4 that the vertical velocity  in blue is close to the log-law when   is big and far when   is small.To quantify this, the errors and relative errors corresponding in the Table 2 show that the bigger   is, the smallest the errors are.
Nevertheless, it is getting steady at some point as we can see in Figure 5, where the errors and the relative errors have been plotted for  ⋆ = 4,  ⋆ = 7 and  ⋆ = 10.The value 10 6 seems to be the threshold value for the coefficient   .It gives for the different values of  ⋆ the curves shown in Figure 6.The peak we can see corresponds to the the height of the viscous sub-layer  0 .We can compare the errors and relative errors between the raw velocity  and the log-law  Log on the one hand, and between the stabilized velocity  + Ψ and  Log on the other hand, for the different regimes given by  ⋆ .We can see in the Table 3 and in the Figure 7 that the stabilization is a better approximation of the log-law, especially around the viscous sub-layer  0 .

With horizontal perturbation
In this part we study the effect of a small oscillation in the horizontal direction, given by () = 0.01 sin   We keep the value   = 10 6 ,  = 1, and we will compare the vertical velocities  we get with the ones we had without this perturbation for the three different  ⋆ .We obtain the Table 4.This shows that even for a small perturbation, the error variation is quite large compared to the values we have.When we add the stabilization function Ψ, we can manage to keep the errors still small, even with the perturbation .This opens an interesting stability problem à la "Lyapunov".

Conclusion and perspectives
We have shown numerically that the solution  of the problem (1.4) is a very good approximation for the known Monin-Obukhov log-law when  = 1 for adiabatic flows, with a large calibration constant   , imposing a Dirichlet condition at Γ top .This is valid for a large range of wind regimes.The numerical code seems to be a good tool for calculating the stabilization function, which could be improved by using a formal mathematical calculation tool, which we have not done yet.
The next step is to couple system (1.4) with the equation for the temperature  , where the source term  is given by (2.7), which is a work under progress.

Figure 2 .
Figure 2. Velocity profiles for different horizontal values.

Figure 5 .
Figure 5. Errors and relative errors dependence on   .

Figure 6 .
Figure 6.Stabilization functions for different values of  ⋆ .Table3.Errors and relative errors between the stabilized/unstabilized velocities and the log-law.
Remark 3.3.Assumption (3.10) about  is not optimal and could be weakened by taking for instance  ∈  −, (Γ top ) for some  > 0,  > 1 depending on , so that it is put in duality with traces on Γ top of test get in this case the existence of a weak solution, we do not know whether the energy inequality (3.13) still holds, or even if it makes sense because of the boundary term ∫︀ .13) Remark 3.2.The existence of a solution still holds when  = 0. Γtop .It seems that there is an interesting theoretical issue at this point.Remark 3.4.The result still holds if one takes the Navier law (2.8) for a given  ∈  2 (Γ top ),  > 0.
Note that the trace of  ∈  1 is not defined at  = 0, where the weight vanishes.On the other hand  |=top is well defined in  1/2 (Γ top ) since the function  belongs to  1,2 in a neighborhood of the top part of the boundary.Hence, the integrals ∫︀ Γtop   and ∫︀ Γtop  are properly defined.

Table 2 .
CD calibration for the different regimes.

Table 4 .
Errors and relative errors between the velocities and the log-law.
⋆ w/o  with  with  and Ψ w/o  with  with  and Ψ