# Surface boundary layers through a scalar equation with an eddy viscosity vanishing at the ground

Luigi C. Berselli, François Legeais, Roger Lewandowski

► To cite this version:

Luigi C. Berselli, François Legeais, Roger Lewandowski. Surface boundary layers through a scalar equation with an eddy viscosity vanishing at the ground. ESAIM: Mathematical Modelling and Numerical Analysis, 2024, 58 (2), pp.489-513. 10.1051/m2an/2024009. hal-04533458v3

HAL Id: hal-04533458

https://hal.science/hal-04533458v3

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Surface boundary layers through a scalar equation with an eddy viscosity vanishing at the ground

Luigi C. Berselli, François Legeais, Roger Lewandowski

► To cite this version:

Luigi C. Berselli, François Legeais, Roger Lewandowski. Surface boundary layers through a scalar equation with an eddy viscosity vanishing at the ground. ESAIM: Mathematical Modelling and Numerical Analysis, 2024, 58 (2), pp.489-513. 10.1051/m2an/2024009. hal-04533458v3

HAL Id: hal-04533458

https://hal.science/hal-04533458v3

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SURFACE BOUNDARY LAYERS THROUGH A SCALAR EQUATION WITH AN EDDY VISCOSITY VANISHING AT THE GROUND

LUIGI C. BERSELLI<fnref n="1" />, FRANÇOIS LEGEAIS<fnref n="2" /> AND ROGER LEWANDOWSKI<fnref n="2" />,*

**Abstract.** We introduce a scalar elliptic equation defined on a boundary layer given by $Π_2 \times [0, z_{top}]$, where $Π_2$ is a two dimensional torus, with an eddy vertical viscosity of order $z^α$, $α \in [0, 1]$, an homogeneous boundary condition at $z = 0$, and a Robin condition at $z = z_{top}$. We show the existence of weak solutions to this boundary problem, distinguishing the cases $0 \le α < 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.

**Mathematics Subject Classification.** 35Q30, 35D30, 76D03, 76D05.

Received January 16, 2023. Accepted January 30, 2024.

## 1. 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 $k-ε$ 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

$$u_{\star} \kappa \frac{d}{dz} \left( z \frac{du}{dz} \right) = 0,$$

together with appropriate boundary conditions, where $u = u(z)$ denotes the mean horizontal velocity component that is supposed to only depend on the height $z > 0$; in this framework, $u_{\star}$ is the friction velocity and $κ$ the

Keywords and phrases. Boundary layers, eddy viscosities, degenerate elliptic equations.

<footnote n="1">¹ Università di Pisa, Dipartimento di Matematica, Via Buonarroti 1/c, I-56127 Pisa, Italy.</footnote>

<footnote n="2">² IRMAR, UMR CNRS 6625, University of Rennes and Odyssey team, INRIA Rennes, France.</footnote>

<footnote>*Corresponding author: roger.lewandowski@univ-rennes1.fr</footnote>
# SURFACE BOUNDARY LAYERS THROUGH A SCALAR EQUATION WITH AN EDDY VISCOSITY VANISHING AT THE GROUND

LUIGI C. BERSELLI<fnref n="1" />, FRANÇOIS LEGEAIS<fnref n="2" /> AND ROGER LEWANDOWSKI<fnref n="2" />,*

**Abstract.** We introduce a scalar elliptic equation defined on a boundary layer given by $Π_2 \times [0, z_{top}]$, where $Π_2$ is a two dimensional torus, with an eddy vertical viscosity of order $z^α$, $α \in [0, 1]$, an homogeneous boundary condition at $z = 0$, and a Robin condition at $z = z_{top}$. We show the existence of weak solutions to this boundary problem, distinguishing the cases $0 \le α < 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.

**Mathematics Subject Classification.** 35Q30, 35D30, 76D03, 76D05.

Received January 16, 2023. Accepted January 30, 2024.

## 1. 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 $k - ε$ 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

$$u_{\star} \kappa \frac{d}{dz} \left( z \frac{du}{dz} \right) = 0,$$

together with appropriate boundary conditions, where $u = u(z)$ denotes the mean horizontal velocity component that is supposed to only depend on the height $z > 0$; in this framework, $u_{\star}$ is the friction velocity and $κ$ the

Keywords and phrases. Boundary layers, eddy viscosities, degenerate elliptic equations.

<footnote n="1">¹ Università di Pisa, Dipartimento di Matematica, Via Buonarroti 1/c, I-56127 Pisa, Italy.</footnote>

<footnote n="2">² IRMAR, UMR CNRS 6625, University of Rennes and Odyssey team, INRIA Rennes, France.</footnote>

<footnote>*Corresponding author: roger.lewandowski@univ-rennes1.fr</footnote>
von Kármán constant. We wonder if a similar simple multi dimensional PDE satisfied by $u = u(\mathbf{x}_h, z)$, $\mathbf{x}_h \in \mathbb{R}^d$ ($d = 1, 2$), is able to yield Monin-Obukhov profiles type, which suggests to introduce the equation,

$$-\nu_h \Delta_h u - \partial_z (\nu_{\text{turb}}(z) \partial_z u) = f, \qquad (1.1)$$

where $f$ is the Boussinesq force specified by a temperature supposed to be constant in this paper, $\nu_{\text{turb}} = \nu_{\text{turb}}(z)$ is an eddy viscosity and $\nu_h > 0$ an horizontal viscosity coefficient. According to standard assumptions and dimensional analysis [5, 10, 15], we should have

$$\nu_{\text{turb}}(z) = \kappa u_{\star} z, \qquad (1.2)$$

in the domain $0 < z < z_{\text{top}}$, where $z_{\text{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–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 $H^{1/2}$ regularity, which does not allow to set a boundary condition at $z = 0$. This is why in many models, $\nu_{\text{turb}}$ is taken to be constant in a viscous sublayer $[0, z_0]$. The same issue occurs in the case of the Smagorinsky model [23], where the eddy viscosity denoted by $\nu_{\text{smag}}$ is given by $\nu_{\text{smag}} = \kappa z^2 |\partial_z u|$ near the ground $z = 0$. This is why in the Smagorinsky case, several authors have suggested to replace the physical $\nu_{\text{smag}}$ by $\nu_{\text{smag}} = \kappa z_0^{2(1-\alpha)} z^{2\alpha} |\partial_z u|$ for some $0 < \alpha < 1$ [4, 20–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

$$\nu_{\text{turb}}(z) = \kappa z_0^{1-\alpha} u_{\star} z^{\alpha}. \qquad (1.3)$$

There is the question of the boundary conditions. It is natural to set $u = 0$ at the ground $z = 0$. Following [9], we take a Navier friction condition at the top of the boundary layer $z = z_{\text{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<fnref n="1" />,

$$\left\{
\begin{array}{ll}
\lambda u - \nu_h \Delta_h u - \mu \partial_z (z^\alpha \partial_z u) = f & \text{in BL,} \\
\mu z_{\text{top}}^\alpha \frac{\partial u(\mathbf{x}_h, z)}{\partial z} = C_D (V(\mathbf{x}_h) - u(\mathbf{x}_h, z_{\text{top}})) & \text{on } \Gamma_{\text{top}}, \\
u = 0 & \text{on } \Gamma_b,
\end{array}
\right.
\qquad (1.4)$$

where $\mathbf{x}_h = (x, y)$ and

$$\text{BL} = \Pi_2 \times [0, z_{\text{top}}], \quad \Gamma_b = \Pi_2 \times \{z=0\}, \quad \Gamma_{\text{top}} = \Pi_2 \times \{z=z_{\text{top}}\}, \qquad (1.5)$$

$0 < \alpha \le 1$, for a given 2D torus denoted by $\Pi_2$, where the term $\lambda u$, $\lambda \ge 0$, is a stabilizing term, useful only in the case $\alpha = 1$, $\mu > 0$ and $C_D > 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 \le \alpha < 1$, see Theorem 3.1 below. Then for $\alpha = 1$, we prove the existence of a very weak solution $u \in H^{1/2}(\text{BL})$ for $\lambda > 0$, in Theorem 4.1, which is the main result of this paper. In this result, we do not impose $u = 0$ at $z = 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 $\alpha = 1$ remains the most accurate to parameterize the SBL, and we are able to calculate

<footnote n="1">¹The acronym SBL has been replaced by BL for simplicity.</footnote>
von Kármán constant. We wonder if a similar simple multi dimensional PDE satisfied by $u = u(\mathbf{x}_h, z)$, $\mathbf{x}_h \in \mathbb{R}^d$ ($d = 1, 2$), is able to yield Monin-Obukhov profiles type, which suggests to introduce the equation,

$$ -\nu_h \Delta_h u - \partial_z (\nu_{\text{turb}}(z) \partial_z u) = f, \qquad (1.1) $$

where $f$ is the Boussinesq force specified by a temperature supposed to be constant in this paper, $\nu_{\text{turb}} = \nu_{\text{turb}}(z)$ is an eddy viscosity and $\nu_h > 0$ an horizontal viscosity coefficient. According to standard assumptions and dimensional analysis [5, 10, 15], we should have

$$ \nu_{\text{turb}}(z) = \kappa u_{\star} z, \qquad (1.2) $$

in the domain $0 < z < z_{\text{top}}$, where $z_{\text{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–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 $H^{1/2}$ regularity, which does not allow to set a boundary condition at $z = 0$. This is why in many models, $\nu_{\text{turb}}$ is taken to be constant in a viscous sublayer $[0, z_0]$. The same issue occurs in the case of the Smagorinsky model [23], where the eddy viscosity denoted by $\nu_{\text{smag}}$ is given by $\nu_{\text{smag}} = \kappa z^2 |\partial_z u|$ near the ground $z = 0$. This is why in the Smagorinsky case, several authors have suggested to replace the physical $\nu_{\text{smag}}$ by $\nu_{\text{smag}} = \kappa z_0^{2(1-\alpha)} z^{2\alpha} |\partial_z u|$ for some $0 < \alpha < 1$ [4, 20–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

$$ \nu_{\text{turb}}(z) = \kappa z_0^{1-\alpha} u_{\star} z^{\alpha}. \qquad (1.3) $$

There is the question of the boundary conditions. It is natural to set $u = 0$ at the ground $z = 0$. Following [9], we take a Navier friction condition at the top of the boundary layer $z = z_{\text{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<fnref n="1" />,

$$ \left\{ \begin{array}{ll} \lambda u - \nu_h \Delta_h u - \mu \partial_z (z^\alpha \partial_z u) = f & \text{in BL,} \\ \mu z_{\text{top}}^\alpha \frac{\partial u(\mathbf{x}_h, z)}{\partial z} = C_D (V(\mathbf{x}_h) - u(\mathbf{x}_h, z_{\text{top}})) & \text{on } \Gamma_{\text{top}}, \\ u = 0 & \text{on } \Gamma_b, \end{array} \right. \qquad (1.4) $$

where $\mathbf{x}_h = (x, y)$ and

$$ \text{BL} = \Pi_2 \times [0, z_{\text{top}}], \quad \Gamma_b = \Pi_2 \times \{z=0\}, \quad \Gamma_{\text{top}} = \Pi_2 \times \{z=z_{\text{top}}\}, \qquad (1.5) $$

$0 < \alpha \le 1$, for a given 2D torus denoted by $\Pi_2$, where the term $\lambda u$, $\lambda \ge 0$, is a stabilizing term, useful only in the case $\alpha = 1$, $\mu > 0$ and $C_D > 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 \le \alpha < 1$, see Theorem 3.1 below. Then for $\alpha = 1$, we prove the existence of a very weak solution $u \in H^{1/2}(\text{BL})$ for $\lambda > 0$, in Theorem 4.1, which is the main result of this paper. In this result, we do not impose $u = 0$ at $z = 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 $\alpha = 1$ remains the most accurate to parameterize the SBL, and we are able to calculate

<footnote n="1">¹The acronym SBL has been replaced by BL for simplicity.</footnote>
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 \le \alpha < 1$, by viscous regularization and proving that the natural weighted space related to the problem is embedded in a standard Sobolev space $W^{1,\gamma}(BL)$. Then we focus on the case $\alpha = 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 $H^{1/2}(BL)$. The last section of the paper is devoted to the numerical results.

## 2. 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

$$ \mathbf{u} = \mathbf{u}(\mathbf{x}_h, z) = (\mathbf{u}_h(\mathbf{x}_h, z), w(\mathbf{x}_h, z)) $$

instead of $\bar{\mathbf{u}}$ for the simplicity, where

$$ \mathbf{u}_h = (u, v), \quad z \in ]0, z_{\text{top}}[$$

$z_{\text{top}} > 0$ being the bottom of (BL). For instance, if (BL) models the surface boundary layer, $z_{\text{top}} \approx 100$ m. We also will need to consider the roughness length $z_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.

### 2.1. Assumptions, general equation and issues

Let $\nu > 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 $w = 0$ and, even if it means making a change of coordinates, we can assume $v = 0$;

- the mean velocity $\mathbf{u} = (u, 0, 0)$ depends only on the altitude, that is $u = u(z)$;

- the eddy viscosity $\nu_{\text{turb}}$ depends on $z$ and $u^* = \sqrt{\nu|\partial_z u(0)|}$, the so-called frictional velocity, which is the tuning parameter of the system (see [18, 24]), which yields

$$ \nu_{\text{turb}} = \nu_{\text{turb}}(z) = C_{\nu}C_{\star}u^{*}z, \qquad (2.1) $$

where $C_{\nu} \approx 15$, $C_{\star} \approx 10$ are non dimensional constants, that we have calibrated by numerical simulations. Typical values of $u_{\star}$ range from 2 to $10 \text{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 $\mathbf{u} = (u, 0, 0)$,

$$ -\frac{d}{dz}\left(\nu_{\text{turb}}(z)\frac{du}{dz}\right) = 0, \qquad (2.2) $$

which formula, once integrated between a given $z_0$ and $z_{\text{top}}$ with appropriate boundary conditions, yields the well known log-law, uniform in $\mathbf{x}_h$, using the calibration constants $C_{\nu}$ and $C_{\star}$:

$$ u(z) = \frac{C_{\star}u^{*}}{C_{\nu}}\left(\log\left(\frac{z}{z_{0}}\right) + 1\right), \quad z \in [z_{0}, z_{\text{top}}]. \qquad (2.3) $$
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 \le \alpha < 1$, by viscous regularization and proving that the natural weighted space related to the problem is embedded in a standard Sobolev space $W^{1,\gamma}(BL)$. Then we focus on the case $\alpha = 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 $H^{1/2}(BL)$. The last section of the paper is devoted to the numerical results.

## 2. 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

$$ \mathbf{u} = \mathbf{u}(\mathbf{x}_h, z) = (\mathbf{u}_h(\mathbf{x}_h, z), w(\mathbf{x}_h, z)) $$

instead of $\bar{\mathbf{u}}$ for the simplicity, where

$$ \mathbf{u}_h = (u, v), \quad z \in ]0, z_{\text{top}}[$$

$z_{\text{top}} > 0$ being the bottom of (BL). For instance, if (BL) models the surface boundary layer, $z_{\text{top}} \approx 100$ m. We also will need to consider the roughness length $z_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.

### 2.1. Assumptions, general equation and issues

Let $\nu > 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 $w = 0$ and, even if it means making a change of coordinates, we can assume $v = 0$;

- the mean velocity $\mathbf{u} = (u, 0, 0)$ depends only on the altitude, that is $u = u(z)$;

- the eddy viscosity $\nu_{\text{turb}}$ depends on $z$ and $u^* = \sqrt{\nu|\partial_z u(0)|}$, the so-called frictional velocity, which is the tuning parameter of the system (see [18, 24]), which yields

$$ \nu_{\text{turb}} = \nu_{\text{turb}}(z) = C_{\nu}C_{\star}u^{*}z, \qquad (2.1) $$

where $C_{\nu} \approx 15$, $C_{\star} \approx 10$ are non dimensional constants, that we have calibrated by numerical simulations. Typical values of $u_{\star}$ range from 2 to $10 \text{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 $\mathbf{u} = (u, 0, 0)$,

$$ -\frac{d}{dz}\left(\nu_{\text{turb}}(z)\frac{du}{dz}\right) = 0, \qquad (2.2) $$

which formula, once integrated between a given $z_0$ and $z_{\text{top}}$ with appropriate boundary conditions, yields the well known log-law, uniform in $\mathbf{x}_h$, using the calibration constants $C_{\nu}$ and $C_{\star}$:

$$ u(z) = \frac{C_{\star}u^{*}}{C_{\nu}}\left(\log\left(\frac{z}{z_{0}}\right) + 1\right), \quad z \in [z_{0}, z_{\text{top}}]. \qquad (2.3) $$
Generally, for $z \in [0, z_0]$, called the viscous sub-layer, a linear profile is considered such that $u = u(z)$ is continuous over $]0, z_{top}]$, and $u(0) = 0$, which means

$$u(z) = \frac{C_{\star} u_{\star}}{C_{\nu}} z, \quad z \in [0, z_0]. \qquad (2.4)$$

Let $u_{\text{Log}}$ denotes the function defined over $]0, z_{\text{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 $u_{\text{Log}}$ to get the right velocity profile, according to the Monin–Obukhov theory [14], which means that

$$u(z) = u_{\text{Log}}(z) + \Psi(\mathbf{x}_h, z), \qquad (2.5)$$

where the function $\Psi(\mathbf{x}_h, z)$ 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, $C_{\nu}$ stands for the von Kármán constant, the value of which being equal to 0.4, and $C_{\star} = 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, z_{\text{top}}[$;

(2) includes an eddy viscosity of the same form as that given by (2.1), where the profile $u = u(\mathbf{x}_h, 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 = $\mathbb{R}^2 \times (0, z_{\text{top}})$:

$$\lambda u - \nu_h \Delta_h u - \partial_z (\nu_{\text{turb}}(z) \partial_z u) = f, \qquad (2.6)$$

for some $\lambda > 0, \nu_h > 0, \Delta_h = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$. The term $\lambda u$ in (2.6) stands for a numerical artefact of an evolutionary term $\partial_t u$, and serves as a system stabilizer, especially in the case $\alpha = 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 $f$ is the Boussinesq force, namely

$$f = g\beta(T_0 - T), \qquad (2.7)$$

where $T$ is the temperature of the fluid, $T_0$ its value at the ground, $g \approx 10 \text{ ms}^{-1}$ is the gravity coefficient, $\beta$ is the coefficient of thermal expansion, a typical value of which for dry air is varies between $0.002 \text{ K}^{-1}$ and $0.003 \text{ K}^{-1}$.

## 2.2. 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 $\Gamma_b$, $u(\mathbf{x}_h, 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

$$C_{\nu} u^{\star} z_{\text{top}} \frac{du}{dz}(\mathbf{x}_{h}, z_{\text{top}}) = C_{D}(V(\mathbf{x}_{h}) - u(\mathbf{x}_{h}, z_{\text{top}})), \qquad (2.8)$$

where $C_D > 0$ is a frictional coefficient and $V = V(\mathbf{x}_h)$ is the velocity of the top layer. In the numerical simulations we have taken $V(\mathbf{x}_h) = u_{\text{Log}}(z_{\text{top}})(1 + \varepsilon(\mathbf{x}_h))$, where $u_{\text{Log}}$ is that given by (2.3) and (2.4). The coefficient $C_D$ will be numerically optimized in function of $u_{\star}, \varepsilon(\mathbf{x}_h)$ is a small perturbation term.

**Remark 2.2.** According to the results of [9], we expect that for large values of $C_D$, the boundary condition (2.8) converges (in some sense) to the continuity condition $u(\mathbf{x}_h) = V(\mathbf{x}_h)$ at $\Gamma_{\text{top}}$, which is well confirmed in this framework by the numerical simulations.
Generally, for $z \in [0, z_0]$, called the viscous sub-layer, a linear profile is considered such that $u = u(z)$ is continuous over $]0, z_{top}]$, and $u(0) = 0$, which means

$$u(z) = \frac{C_{\star} u_{\star}}{C_{\nu}} z, \quad z \in [0, z_0]. \qquad (2.4)$$

Let $u_{\text{Log}}$ denotes the function defined over $]0, z_{\text{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 $u_{\text{Log}}$ to get the right velocity profile, according to the Monin–Obukhov theory [14], which means that

$$u(z) = u_{\text{Log}}(z) + \Psi(\mathbf{x}_h, z), \qquad (2.5)$$

where the function $\Psi(\mathbf{x}_h, z)$ 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, $C_{\nu}$ stands for the von Kármán constant, the value of which being equal to 0.4, and $C_{\star} = 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, z_{\text{top}}[$;

(2) includes an eddy viscosity of the same form as that given by (2.1), where the profile $u = u(\mathbf{x}_h, 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 = $\mathbb{R}^2 \times (0, z_{\text{top}})$:

$$\lambda u - \nu_h \Delta_h u - \partial_z (\nu_{\text{turb}}(z) \partial_z u) = f, \qquad (2.6)$$

for some $\lambda > 0, \nu_h > 0, \Delta_h = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$. The term $\lambda u$ in (2.6) stands for a numerical artefact of an evolutionary term $\partial_t u$, and serves as a system stabilizer, especially in the case $\alpha = 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 $f$ is the Boussinesq force, namely

$$f = g\beta(T_0 - T), \qquad (2.7)$$

where $T$ is the temperature of the fluid, $T_0$ its value at the ground, $g \approx 10 \text{ ms}^{-1}$ is the gravity coefficient, $\beta$ is the coefficient of thermal expansion, a typical value of which for dry air is varies between $0.002 \text{ K}^{-1}$ and $0.003 \text{ K}^{-1}$.

## 2.2. 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 $\Gamma_b$, $u(\mathbf{x}_h, 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

$$C_{\nu} u^{*} z_{\text{top}} \frac{du}{dz}(\mathbf{x}_{h}, z_{\text{top}}) = C_{D}(V(\mathbf{x}_{h}) - u(\mathbf{x}_{h}, z_{\text{top}})), \qquad (2.8)$$

where $C_D > 0$ is a frictional coefficient and $V = V(\mathbf{x}_h)$ is the velocity of the top layer. In the numerical simulations we have taken $V(\mathbf{x}_h) = u_{\text{Log}}(z_{\text{top}})(1 + \varepsilon(\mathbf{x}_h))$, where $u_{\text{Log}}$ is that given by (2.3) and (2.4). The coefficient $C_D$ will be numerically optimized in function of $u_\star$, $\varepsilon(\mathbf{x}_h)$ is a small perturbation term.

**Remark 2.2.** According to the results of [9], we expect that for large values of $C_D$, the boundary condition (2.8) converges (in some sense) to the continuity condition $u(\mathbf{x}_h) = V(\mathbf{x}_h)$ at $\Gamma_{\text{top}}$, which is well confirmed in this framework by the numerical simulations.
## 2.3. 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 $H^{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 $\nu_{\text{turb}}$ of the form

$$ \nu_{\text{turb}} = \nu_{\text{turb}}(z) = z_0^{1-\alpha} C_{\nu} u^* z^{\alpha}, \quad \alpha \in [0, 1], \qquad (2.9) $$

the main feature of which is that it degenerates at the ground, but with a different velocity. The parameter $z_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:

$$ \left\{ \begin{array}{l} C_D u(\mathbf{x}_h, z_{\text{top}}) + \mu z_{\text{top}}^{\alpha} \frac{\partial u}{\partial z}(\mathbf{x}_h, z_{\text{top}}) = G(\mathbf{x}_h) \\ u(x_h, 0) = 0, \end{array} \right. \qquad (2.10) $$

for $\mu = z_0^{1-\alpha} C_{\nu} u^*$ and where $G = C_D V$.

It remains to clarify the boundary conditions in the $\mathbf{x}_h$-axis. For practical numerical simulations, we have to limit ourselves to a finite computational box $[0, L_x] \times [0, L_y] \times 0$, $z_{\text{top}}$ [for given scales $L_x$ and $L_y$, which raises the question of the boundary conditions at the entrance, exit and sides of the computational box, namely

$$ \begin{aligned} \Gamma_{\text{in}} &= \{0\} \times [0, L_y] \times [0, z_{\text{top}}], & \Gamma_{\text{out}} &= \{L_x\} \times [0, L_y] \times [0, z_{\text{top}}], \\ \Gamma_{s,1} &= [0, L_x] \times \{0\} \times 0, & \Gamma_{s,2} &= [0, L_x] \times \{L_y\} \times [0, z_{\text{top}}]. \end{aligned} $$

In [11] we have considered at $\Gamma_{\text{in}}$ a fixed given field in $H_{00}^{1/2}(\Gamma_{\text{in}})$ and nonlinear Neumann transparent boundary conditions at $\Gamma_0$ which yields serious technical issues, both theoretically and numerically. In this paper, we opt for horizontal periodic boundary conditions, which means that $u$ must satisfy for all formal derivative operator $D^n$, $n \ge 0$,

$$ D^n u(x + L_x, y + L_y, z) = D^n u(x, y, z), \quad \forall \mathbf{x}_h = (x, y) \in \mathbb{R}^2, \forall z \in [0, z_{\text{top}}], \qquad (2.11) $$

and imposing the invariance for translations implies working with the torus $\Pi_2$ given by

$$ \Pi_2 = ([0, L_x] \times [0, L_y]) / \mathcal{T}_2, \qquad (2.12) $$

where $\mathcal{T}_2$ denotes the set of wave vectors given by:

$$ \mathcal{T}_2 = \frac{2\pi}{L_x} \mathbb{Z} \times \frac{2\pi}{L_y} \mathbb{Z}. \qquad (2.13) $$

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.

## 3. THE $0 < \alpha < 1$ CASE

Throughout this section, we assume that $0 < \alpha < 1$ and $\lambda > 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.
## 2.3. 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 $H^{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 $\nu_{\text{turb}}$ of the form

$$ \nu_{\text{turb}} = \nu_{\text{turb}}(z) = z_0^{1-\alpha} C_{\nu} u^* z^{\alpha}, \quad \alpha \in [0, 1], \qquad (2.9) $$

the main feature of which is that it degenerates at the ground, but with a different velocity. The parameter $z_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:

$$ \left\{ \begin{array}{l} C_D u(\mathbf{x}_h, z_{\text{top}}) + \mu z_{\text{top}}^{\alpha} \frac{\partial u}{\partial z}(\mathbf{x}_h, z_{\text{top}}) = G(\mathbf{x}_h) \\ u(x_h, 0) = 0, \end{array} \right. \qquad (2.10) $$

for $\mu = z_0^{1-\alpha} C_{\nu} u^*$ and where $G = C_D V$.

It remains to clarify the boundary conditions in the $\mathbf{x}_h$-axis. For practical numerical simulations, we have to limit ourselves to a finite computational box $[0, L_x] \times [0, L_y] \times 0$, $z_{\text{top}}$ [for given scales $L_x$ and $L_y$, which raises the question of the boundary conditions at the entrance, exit and sides of the computational box, namely

$$ \begin{array}{ll} \Gamma_{\text{in}} = \{0\} \times [0, L_y] \times [0, z_{\text{top}}], & \Gamma_{\text{out}} = \{L_x\} \times [0, L_y] \times [0, z_{\text{top}}], \\ \Gamma_{s,1} = [0, L_x] \times \{0\} \times 0, & \Gamma_{s,2} = [0, L_x] \times \{L_y\} \times [0, z_{\text{top}}]. \end{array} $$

In [11] we have considered at $\Gamma_{\text{in}}$ a fixed given field in $H_{00}^{1/2}(\Gamma_{\text{in}})$ and nonlinear Neumann transparent boundary conditions at $\Gamma_0$ which yields serious technical issues, both theoretically and numerically. In this paper, we opt for horizontal periodic boundary conditions, which means that $u$ must satisfy for all formal derivative operator $D^n$, $n \ge 0$,

$$ D^n u(x + L_x, y + L_y, z) = D^n u(x, y, z), \quad \forall \mathbf{x}_h = (x, y) \in \mathbb{R}^2, \forall z \in [0, z_{\text{top}}], \qquad (2.11) $$

and imposing the invariance for translations implies working with the torus $\Pi_2$ given by

$$ \Pi_2 = ([0, L_x] \times [0, L_y]) / \mathcal{T}_2, \qquad (2.12) $$

where $\mathcal{T}_2$ denotes the set of wave vectors given by:

$$ \mathcal{T}_2 = \frac{2\pi}{L_x} \mathbb{Z} \times \frac{2\pi}{L_y} \mathbb{Z}. \qquad (2.13) $$

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.

## 3. THE $0 < \alpha < 1$ CASE

Throughout this section, we assume that $0 < \alpha < 1$ and $\lambda > 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

$$ BL = \Pi_2 \times [0, z_{\text{top}}], \quad \Gamma_b = \Pi_2 \times \{z=0\}, \quad \Gamma_{\text{top}} = \Pi_2 \times \{z=z_{\text{top}}\}. \qquad (3.1) $$

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:

$$ \lambda \int_{BL} u^2 + \nu_h \int_{BL} |\nabla_h u|^2 + \mu \int_{BL} z^\alpha (\partial_z u)^2 + C_D \int_{\Gamma_{\text{top}}} |u|^2 - \int_{\Gamma_{\text{top}}} uG = \int_{BL} fu. \qquad (3.2) $$

## 3.1. 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 $D_b(\text{BL})$ denote the space of functions $u \in C^\infty(\text{BL})$ such that there exists $\delta = \delta(u) \in ]0, 1[$, ensuring that the support of $u$ is contained within $\Pi_2 \times [\delta(u), 1]$. This space is endowed with the standard Schwartz topology, characterized by the semi-norms $\|D^n u\|_\infty$. In simpler terms, $D_b(\text{BL})$ consists of $C^\infty$ functions that vanish in a neighborhood of the bottom $\Gamma_b$ of the boundary layer.

- Let $D'_b(\text{BL})$ denote the topological dual of $D_b(\text{BL})$, a sort of distributional space that inspires our chosen notation.

- $W_{0,b}^{1,\gamma}(\text{BL})$ ($\gamma > 1$) denotes the closure of $D_b(\text{BL})$ for the norm $\|\nabla_h u\|_{0,\gamma} + \|\partial_z u\|_{0,\gamma}$, where $\|\cdot\|_{s,p}$ stands for the usual $W^{s,p}$ norm.

In particular,

$$ W_{0,b}^{1,\gamma}(\text{BL}) = \{u \in W^{1,\gamma}(\text{BL}), u = 0 \text{ on } \Gamma_b\}. \qquad (3.3) $$

According to the energy balance (3.2), we are led to consider the natural weighted space $V_\alpha$ defined as the closure of $D_b(\text{BL})$ equipped with the norm

$$ \forall u \in D_b(\text{BL}), \quad \|u\|_{V_\alpha} = \left( \|\nabla_h u\|_{0,2}^2 + \|z^{\alpha/2} \partial_z u\|_{0,2}^2 \right)^{1/2}. \qquad (3.4) $$

**Proposition 3.1.** For all $\alpha \in [0, 1[$ we have the embedding,

$$ \forall \gamma \in ]1, \frac{2}{\alpha+1}] \quad , \quad V_\alpha \subset W_{0,b}^{1,\gamma}(\text{BL}), $$

and the following inequality holds $\forall u \in V_\alpha$,

$$ \|u\|_{1,\gamma} \le C \left( \frac{2-\gamma}{2-(1+\alpha)\gamma} \right)^{\frac{2}{\gamma}-1} \|u\|_{V_\alpha}, \qquad (3.5) $$

for some constant $C > 0$.

*Proof.* Let $\gamma > 0$ and $\rho > 0$ that will be fixed later. Let $u \in D_{0,b}(\text{BL})$. The Hölder inequality yields:

$$ \int_0^1 |\partial_z u|^\gamma dz = \int_0^1 z^\rho |\partial_z u|^\gamma \frac{1}{z^\rho} dz \le \left( \int_0^1 z^{2\rho/\gamma} |\partial_z u|^2 dz \right)^{\gamma/2} \left( \int_0^1 \frac{dz}{z^{2\rho/2-\gamma}} \right)^{\frac{2-\gamma}{2}}, $$

with $\frac{\gamma}{2} + \frac{2-\gamma}{2} = 1$. The second integral is well-defined if and only if $\frac{2\rho}{2-\gamma} < 1$. Then, choosing $\rho$ such that $\frac{2\rho}{\gamma} = \alpha$ yields the condition $\gamma < \frac{2}{\alpha+1}$. Inequality (3.5) follows after an elementary calculation and integration with respect to the $dx_h$ variables. $\square$

It follows from Proposition 3.1 and standard reasoning on Sobolev spaces, that functions in $V_\alpha$ have a trace at $z = 0$ equal to zero, and also we have the following characterization

$$ V_\alpha = \{u \in D'_b(\text{BL}) \text{ s.t. } \nabla_h u \in L^2(\text{BL}), z^{\alpha/2} \partial_z u \in L^2(\text{BL}), u = 0 \text{ on } \Gamma_b\}. \qquad (3.6) $$
We will use standard Lebesgue and Sobolev spaces and from now on

$$ BL = \Pi_2 \times [0, z_{\text{top}}], \quad \Gamma_b = \Pi_2 \times \{z=0\}, \quad \Gamma_{\text{top}} = \Pi_2 \times \{z=z_{\text{top}}\}. \qquad (3.1) $$

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:

$$ \lambda \int_{BL} u^2 + \nu_h \int_{BL} |\nabla_h u|^2 + \mu \int_{BL} z^{\alpha} (\partial_z u)^2 + C_D \int_{\Gamma_{\text{top}}} |u|^2 - \int_{\Gamma_{\text{top}}} uG = \int_{BL} fu. \qquad (3.2) $$

## 3.1. 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 $D_b(BL)$ denote the space of functions $u \in C^\infty(BL)$ such that there exists $\delta = \delta(u) \in ]0, 1[$, ensuring that the support of $u$ is contained within $\Pi_2 \times [\delta(u), 1]$. This space is endowed with the standard Schwartz topology, characterized by the semi-norms $\|D^n u\|_\infty$. In simpler terms, $D_b(BL)$ consists of $C^\infty$ functions that vanish in a neighborhood of the bottom $\Gamma_b$ of the boundary layer.

- Let $D'_b(BL)$ denote the topological dual of $D_b(BL)$, a sort of distributional space that inspires our chosen notation.

- $W_{0,b}^{1,\gamma}(BL)$ ($\gamma > 1$) denotes the closure of $D_b(BL)$ for the norm $\|\nabla_h u\|_{0,\gamma} + \|\partial_z u\|_{0,\gamma}$, where $\|\cdot\|_{s,p}$ stands for the usual $W^{s,p}$ norm.

In particular,

$$ W_{0,b}^{1,\gamma}(BL) = \{u \in W^{1,\gamma}(BL), u = 0 \text{ on } \Gamma_b\}. \qquad (3.3) $$

According to the energy balance (3.2), we are led to consider the natural weighted space $V_\alpha$ defined as the closure of $D_b(BL)$ equipped with the norm

$$ \forall u \in D_b(BL), \quad \|u\|_{V_\alpha} = \left( \|\nabla_h u\|_{0,2}^2 + \|z^{\alpha/2} \partial_z u\|_{0,2}^2 \right)^{1/2}. \qquad (3.4) $$

**Proposition 3.1.** For all $\alpha \in [0, 1[$ we have the embedding,

$$ \forall \gamma \in ]1, \frac{2}{\alpha+1}], \quad V_\alpha \subset W_{0,b}^{1,\gamma}(BL), $$

and the following inequality holds $\forall u \in V_\alpha$,

$$ \|u\|_{1,\gamma} \le C \left( \frac{2-\gamma}{2-(1+\alpha)\gamma} \right)^{\frac{2}{\gamma}-1} \|u\|_{V_\alpha}, \qquad (3.5) $$

for some constant $C > 0$.

**Proof.** Let $\gamma > 0$ and $\rho > 0$ that will be fixed later. Let $u \in D_{0,b}(BL)$. The Hölder inequality yields:

$$ \int_0^1 |\partial_z u|^\gamma dz = \int_0^1 z^\rho |\partial_z u|^\gamma \frac{1}{z^\rho} dz \le \left( \int_0^1 z^{2\rho/\gamma} |\partial_z u|^2 dz \right)^{\gamma/2} \left( \int_0^1 \frac{dz}{z^{2\rho/2-\gamma}} \right)^{\frac{2-\gamma}{2}}, $$

with $\frac{\gamma}{2} + \frac{2-\gamma}{2} = 1$. The second integral is well-defined if and only if $\frac{2\rho}{2-\gamma} < 1$. Then, choosing $\rho$ such that $\frac{2\rho}{\gamma} = \alpha$ yields the condition $\gamma < \frac{2}{\alpha+1}$. Inequality (3.5) follows after an elementary calculation and integration with respect to the $dx_h$ variables. $\square$

It follows from Proposition 3.1 and standard reasoning on Sobolev spaces, that functions in $V_\alpha$ have a trace at $z = 0$ equal to zero, and also we have the following characterization

$$ V_\alpha = \{u \in D'_b(BL) \text{ s.t. } \nabla_h u \in L^2(BL), z^{\alpha/2} \partial_z u \in L^2(BL), u = 0 \text{ on } \Gamma_b\}. \qquad (3.6) $$
## 3.2. Weak formulation

Proposition 3.1 can be rephrased, to work in standard (unweighted) Sobolev spaces, as follows

$$V_{\alpha} \subset W_{0,b}^{1,(\frac{2}{1+\alpha})^{-}} := \bigcap_{1<\gamma<\frac{2}{1+\alpha}} W_{0,b}^{1,\gamma}(\text{BL}), \quad (3.7)$$

which is put in duality with the set

$$W_{0,b}^{1,(\frac{2}{1-\alpha})^{+}} = \bigcup_{\eta > \frac{2}{1-\alpha}} W_{0,b}^{1,\eta}(\text{BL}). \quad (3.8)$$

Throughout the section, we assume that

$$f \in L^2(\text{BL}), \quad (3.9)$$

$$G \in L^2(\Gamma_{\text{top}}). \quad (3.10)$$

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.

**Definition 3.1.** For $\alpha \in [0, 1[$ we say that $u \in V_{\alpha}$ is a weak solution to problem (1.4), if $\forall v \in W_{0}^{1,(\frac{2}{1-\alpha})^{+}}$,

$$\lambda \int_{\text{BL}} uv + \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z^{\alpha} \partial_z u \partial_z v + \int_{\Gamma_{\text{top}}} (C_D u - G)v = \int_{\text{BL}} fv. \quad (3.11)$$

**Remark 3.1.** Note that all the terms in the integrals written in (3.11) are well-defined. However, the solution $u$ 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 < \gamma < \frac{2}{1+\alpha}$ the following estimate in $W_{0,b}^{1,\gamma}(\text{BL})$, satisfied by any given regular solution $u$ to the variational problem (3.11):

$$\|u\|_{0,1,\gamma} \le \frac{C_{\gamma}}{\inf\{\nu_h, \mu\}} (\|f\|_{0,2} + \|G\|_{1-\frac{1}{\gamma},\gamma;\Gamma_{\text{top}}}), \quad (3.12)$$

where $C_{\gamma} \to \infty$ as $\gamma \to \frac{2}{1+\alpha}$.

The aim of the rest of this section is proving the following existence result.

**Theorem 3.1.** Problem (1.4) admits a weak solution $u \in V_{\alpha}$. Moreover, the solution satisfies the energy inequality

$$\lambda \int_{\text{BL}} u^2 + \nu_h \int_{\text{BL}} |\nabla_h u|^2 + \mu \int_{\text{BL}} z^{\alpha} |\partial_z u|^2 + C_D \int_{\Gamma_{\text{top}}} |u|^2 - \int_{\Gamma_{\text{top}}} G u \le \int_{\text{BL}} fu. \quad (3.13)$$

**Remark 3.2.** The existence of a solution still holds when $\lambda = 0$.

**Remark 3.3.** Assumption (3.10) about $G$ is not optimal and could be weakened by taking for instance $G \in W^{-s,p}(\Gamma_{\text{top}})$ for some $s > 0$, $p > 1$ depending on $\alpha$, so that it is put in duality with traces on $\Gamma_{\text{top}}$ of test functions in $W_{0,b}^{1,(\frac{2}{1-\alpha})^{+}}$.

If we still 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 $\int_{\Gamma_{\text{top}}} uG$. 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 $V \in L^2(\Gamma_{\text{top}})$, $\alpha > 0$.
## 3.2. Weak formulation

Proposition 3.1 can be rephrased, to work in standard (unweighted) Sobolev spaces, as follows

$$V_{\alpha} \subset W_{0,b}^{1,(\frac{2}{1+\alpha})^{-}} := \bigcap_{1<\gamma<\frac{2}{1+\alpha}} W_{0,b}^{1,\gamma}(\text{BL}), \qquad (3.7)$$

which is put in duality with the set

$$W_{0,b}^{1,(\frac{2}{1-\alpha})^{+}} = \bigcup_{\eta > \frac{2}{1-\alpha}} W_{0,b}^{1,\eta}(\text{BL}). \qquad (3.8)$$

Throughout the section, we assume that

$$f \in L^2(\text{BL}), \qquad (3.9)$$

$$G \in L^2(\Gamma_{\text{top}}). \qquad (3.10)$$

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.

**Definition 3.1.** For $\alpha \in [0, 1[$ we say that $u \in V_{\alpha}$ is a weak solution to problem (1.4), if $\forall v \in W_{0}^{1,(\frac{2}{1-\alpha})^{+}}$,

$$\lambda \int_{\text{BL}} uv + \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z^{\alpha} \partial_z u \partial_z v + \int_{\Gamma_{\text{top}}} (C_D u - G)v = \int_{\text{BL}} fv. \qquad (3.11)$$

**Remark 3.1.** Note that all the terms in the integrals written in (3.11) are well-defined. However, the solution $u$ 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 < \gamma < \frac{2}{1+\alpha}$ the following estimate in $W_{0,b}^{1,\gamma}(\text{BL})$, satisfied by any given regular solution $u$ to the variational problem (3.11):

$$\|u\|_{0,1,\gamma} \le \frac{C_{\gamma}}{\inf\{\nu_h, \mu\}} (\|f\|_{0,2} + \|G\|_{1-\frac{1}{\gamma},\gamma;\Gamma_{\text{top}}}), \qquad (3.12)$$

where $C_{\gamma} \to \infty$ as $\gamma \to \frac{2}{1+\alpha}$.

The aim of the rest of this section is proving the following existence result.

**Theorem 3.1.** Problem (1.4) admits a weak solution $u \in V_{\alpha}$. Moreover, the solution satisfies the energy inequality

$$\lambda \int_{\text{BL}} u^2 + \nu_h \int_{\text{BL}} |\nabla_h u|^2 + \mu \int_{\text{BL}} z^{\alpha} |\partial_z u|^2 + C_D \int_{\Gamma_{\text{top}}} |u|^2 - \int_{\Gamma_{\text{top}}} G u \le \int_{\text{BL}} fu. \qquad (3.13)$$

**Remark 3.2.** The existence of a solution still holds when $\lambda = 0$.

**Remark 3.3.** Assumption (3.10) about $G$ is not optimal and could be weakened by taking for instance $G \in W^{-s,p}(\Gamma_{\text{top}})$ for some $s > 0$, $p > 1$ depending on $\alpha$, so that it is put in duality with traces on $\Gamma_{\text{top}}$ of test functions in $W_{0,b}^{1,(\frac{2}{1-\alpha})^{+}}$.

If we still 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 $\int_{\Gamma_{\text{top}}} uG$. 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 $V \in L^2(\Gamma_{\text{top}})$, $\alpha > 0$.
### 3.3. Viscous regularization and proof of the existence result

From now we set

$$ \gamma^* = \frac{2}{1+\alpha}, \qquad (3.14) $$

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 z-direction, which means that we consider the following problem, for a given $\varepsilon > 0$:

$$ \left\{ \begin{array}{ll} \lambda u - \nu_h \Delta_h u - \mu \partial_z (z^\alpha \partial_z u) - \epsilon \partial_{zz}^2 u = f, & \text{on } \Gamma_{\text{top}}, \\ C_D u + \mu z_{\text{top}}^\alpha \frac{\partial u}{\partial z} = G & \text{on } \Gamma_{\text{b}}, \\ u = 0 & \text{on } \Gamma_{\text{b}}. \end{array} \right. \qquad (3.15) $$

**Definition 3.2.** Let $\varepsilon > 0$. We say that $u_\varepsilon \in W_{0,b}^{1,2}(\text{BL})$ is a weak solution to problem (3.15), if $\forall v \in W_{0,b}^{1,2}(\text{BL})$,

$$ \lambda \int_{\text{BL}} u_{\varepsilon} v + \nu_h \int_{\text{BL}} \nabla_h u_{\varepsilon} \cdot \nabla_h v + \mu \int_{\text{BL}} z^{\alpha} \partial_z u_{\varepsilon} \partial_z v + \varepsilon \int_{\text{BL}} \partial_z u_{\varepsilon} \partial_z v + C_D \int_{\Gamma_{\text{top}}} u_{\varepsilon} v - \int_{\Gamma_{\text{top}}} G v = \int_{\text{BL}} f v. \qquad (3.16) $$

The existence and uniqueness of a weak solution to problem (3.16) is straightforward by the Lax-Milgram theorem. Moreover, as $u_\varepsilon$ can be taken as test function, it satisfies the following energy balance (equality):

$$ \underbrace{\lambda \int_{\text{BL}} u_{\varepsilon}^2}_{I_{1,\varepsilon}} + \underbrace{\nu_h \int_{\text{BL}} |\nabla_h u_{\varepsilon}|^2 + \int_{\text{BL}} (\varepsilon + \mu z^{\alpha}) |\partial_z u_{\varepsilon}|^2}_{I_{2,\varepsilon}} + \underbrace{C_D \int_{\Gamma_{\text{top}}} |u_{\varepsilon}|^2}_{I_{3,\varepsilon}} - \underbrace{\int_{\Gamma_{\text{top}}} G u_{\varepsilon}}_{I_{4,\varepsilon}} = \underbrace{\int_{\text{BL}} f u_{\varepsilon}}_{I_{5,\varepsilon}}. \qquad (3.17) $$

From this, we are able to finish the proof of Theorem 3.1 by taking the limit in (3.16) when $\varepsilon \to 0$.

We deduce from (3.17) and standard calculus inequalities that the family $(u_\varepsilon)_{\varepsilon>0}$ is uniformly bounded in $V_\alpha$, as well as in $W_{0,b}^{1,\gamma}(\text{BL})$ for any $1 < \gamma < \gamma^*$ say

$$ \|u_{\varepsilon}\|_{0;1,\gamma} \le C(\gamma, G, f, \nu_h, \mu). \qquad (3.18) $$

Arguing as in Chapter 7 of [24], we can extract a (sub)sequence $(u_{\varepsilon_n})_{n \in \mathbb{N}}$ that weakly converges to some $u \in V_\alpha$, which is also weakly converging in $W_{0,b}^{1,\gamma}(\text{BL})$ for all $1 < \gamma < \gamma^*$, and which is strongly converging in $L^2(\text{BL})$.

Moreover, by the trace theorem and the Sobolev theorem, the sequence $(\varepsilon_n)$ can be chosen such that in addition $(\text{tr}[u_{\varepsilon_n}])_{n \in \mathbb{N}}$ strongly converges to $\text{tr}[u]$ in $L^2(\Gamma_{\text{top}})$. Let $1 < \gamma < \gamma^*$ and take as test function

$$ v \in W_{0,b}^{1,\gamma'}(\text{BL}) \subset W_{0,b}^{1,\frac{2}{1-\alpha}} \subset W_{0,b}^{1,2}(\text{BL}), $$

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 $w \in V_\alpha$, and let the linear form $\Psi_v$ given by

$$ \langle \Psi_v, w \rangle = \nu_h \int_{\text{BL}} \nabla_h w \cdot \nabla_h v \, dx + \mu \int_{\text{BL}} z^{\alpha} \partial_z w \partial_z v \, dx. $$

By the Hölder inequality we obtain

$$ |\langle \Psi_v, w \rangle| \le \sup\{\nu_h, \mu z_{\text{top}}^{\alpha}\} \|w\|_{0;1,\gamma} \|v\|_{0;1,\gamma'}, $$
### 3.3. Viscous regularization and proof of the existence result

From now we set

$$ \gamma^* = \frac{2}{1+\alpha}, \qquad (3.14) $$

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 z-direction, which means that we consider the following problem, for a given $\varepsilon > 0$:

$$ \left\{ \begin{array}{ll} \lambda u - \nu_h \Delta_h u - \mu \partial_z (z^\alpha \partial_z u) - \epsilon \partial_{zz}^2 u = f, & \text{on } \Gamma_{\text{top}}, \\ C_D u + \mu z_{\text{top}}^\alpha \frac{\partial u}{\partial z} = G & \text{on } \Gamma_{\text{b}}, \\ u = 0 & \text{on } \Gamma_{\text{b}}. \end{array} \right. \qquad (3.15) $$

**Definition 3.2.** Let $\varepsilon > 0$. We say that $u_\varepsilon \in W_{0,b}^{1,2}(\text{BL})$ is a weak solution to problem (3.15), if $\forall v \in W_{0,b}^{1,2}(\text{BL})$,

$$ \lambda \int_{\text{BL}} u_{\varepsilon} v + \nu_h \int_{\text{BL}} \nabla_h u_{\varepsilon} \cdot \nabla_h v + \mu \int_{\text{BL}} z^{\alpha} \partial_z u_{\varepsilon} \partial_z v + \varepsilon \int_{\text{BL}} \partial_z u_{\varepsilon} \partial_z v + C_D \int_{\Gamma_{\text{top}}} u_{\varepsilon} v - \int_{\Gamma_{\text{top}}} G v = \int_{\text{BL}} f v. \qquad (3.16) $$

The existence and uniqueness of a weak solution to problem (3.16) is straightforward by the Lax-Milgram theorem. Moreover, as $u_\varepsilon$ can be taken as test function, it satisfies the following energy balance (equality):

$$ \underbrace{\lambda \int_{\text{BL}} u_{\varepsilon}^2}_{I_{1,\varepsilon}} + \underbrace{\nu_h \int_{\text{BL}} |\nabla_h u_{\varepsilon}|^2 + \int_{\text{BL}} (\varepsilon + \mu z^{\alpha}) |\partial_z u_{\varepsilon}|^2}_{I_{2,\varepsilon}} + \underbrace{C_D \int_{\Gamma_{\text{top}}} |u_{\varepsilon}|^2}_{I_{3,\varepsilon}} - \underbrace{\int_{\Gamma_{\text{top}}} G u_{\varepsilon}}_{I_{4,\varepsilon}} = \underbrace{\int_{\text{BL}} f u_{\varepsilon}}_{I_{5,\varepsilon}}. \qquad (3.17) $$

From this, we are able to finish the proof of Theorem 3.1 by taking the limit in (3.16) when $\varepsilon \to 0$.

We deduce from (3.17) and standard calculus inequalities that the family $(u_\varepsilon)_{\varepsilon>0}$ is uniformly bounded in $V_\alpha$, as well as in $W_{0,b}^{1,\gamma}(\text{BL})$ for any $1 < \gamma < \gamma^*$ say

$$ \|u_\varepsilon\|_{0;1,\gamma} \le C(\gamma, G, f, \nu_h, \mu). \qquad (3.18) $$

Arguing as in Chapter 7 of [24], we can extract a (sub)sequence $(u_{\varepsilon_n})_{n \in \mathbb{N}}$ that weakly converges to some $u \in V_\alpha$, which is also weakly converging in $W_{0,b}^{1,\gamma}(\text{BL})$ for all $1 < \gamma < \gamma^*$, and which is strongly converging in $L^2(\text{BL})$.

Moreover, by the trace theorem and the Sobolev theorem, the sequence $(\varepsilon_n)$ can be chosen such that in addition $(\text{tr}[u_{\varepsilon_n}])_{n \in \mathbb{N}}$ strongly converges to $\text{tr}[u]$ in $L^2(\Gamma_{\text{top}})$. Let $1 < \gamma < \gamma^*$ and take as test function

$$ v \in W_{0,b}^{1,\gamma'}(\text{BL}) \subset W_{0,b}^{1,\frac{2}{1-\alpha}} \subset W_{0,b}^{1,2}(\text{BL}), $$

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 $w \in V_\alpha$, and let the linear form $\Psi_v$ given by

$$ \langle \Psi_v, w \rangle = \nu_h \int_{\text{BL}} \nabla_h w \cdot \nabla_h v \, dx + \mu \int_{\text{BL}} z^{\alpha} \partial_z w \partial_z v \, dx. $$

By the Hölder inequality we obtain

$$ |\langle \Psi_v, w \rangle| \le \sup\{\nu_h, \mu z_{\text{top}}^{\alpha}\} \|w\|_{0;1,\gamma} \|v\|_{0;1,\gamma'}, $$
therefore $\Psi_v \in W_{0,b}^{1,\gamma}(\text{BL})'$, which leads to

$$ \lim_{n \to \infty} \langle \Psi_v, u_{\epsilon_n} \rangle = \langle \Psi_v, u \rangle = \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z^\alpha \partial_z u \partial_z v. $$

Moreover,

$$ \left| \epsilon_n \int_{\text{BL}} \partial_z u_{\epsilon_n} \partial_z v \right| \le \epsilon_n \|u_{\epsilon_n}\|_{0;1,\gamma} \|v\|_{0;1,\gamma'} \le C(\gamma, G, f, \nu_h, \mu) \epsilon_n \|v\|_{0;1,\gamma'}, $$

giving

$$ \lim_{n \to \infty} \epsilon_n \int_{\text{BL}} \partial_z u_{\epsilon_n} \partial_z v = 0. $$

In addition, considering the properties of the sequence $(u_{\epsilon_n})_{n \in \mathbb{N}}$, we have

$$ \lim_{n \to \infty} \lambda \int_{\text{BL}} u_{\epsilon_n} v = \lambda \int_{\text{BL}} uv, \quad \lim_{n \to \infty} C_D \int_{\Gamma_{\text{top}}} u_{\epsilon_n} v = C_D \int_{\Gamma_{\text{top}}} uv. $$

Therefore $u$ satisfies (3.11).

It remains to show that the energy inequality (3.13) holds. Starting from (3.17), we have on one hand

$$ \lim_{n \to \infty} I_{1, \epsilon_n} = \lambda \int_{\text{BL}} |u|^2, \quad \lim_{n \to \infty} I_{3, \epsilon_n} = C_D \int_{\Gamma_{\text{top}}} |u|^2, $$

$$ \lim_{n \to \infty} I_{4, \epsilon_n} = \int_{\Gamma_{\text{top}}} u G, \quad \lim_{n \to \infty} I_{5, \epsilon_n} = \int_{\text{BL}} f u, $$

and on the other hand by lower semi-continuity

$$ \nu_h \int_{\text{BL}} |\nabla_h u|^2 + \mu \int_{\text{BL}} z^\alpha |\partial_z u|^2 \le \liminf_{n \to \infty} I_{2, \epsilon_n}, $$

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 $V_\alpha$, 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.

### 4. THE CASE $\alpha = 1$

We now consider the system:

$$ \begin{cases} \lambda u - \nu_h \Delta_h u - \mu \partial_z (z \partial_z u) = f & \text{in BL,} \\ C_D u + \mu z_{\text{top}} \frac{\partial u}{\partial z} = G & \text{on } \Gamma_{\text{top}}. \end{cases} \quad (4.1) $$

As we shall see, the best we can get is an estimate in an $H^{1/2}$-like space near the bottom $\Gamma_b$. Therefore, we are not able to define the trace of such functions at $\Gamma_b$. This is why we do not impose any boundary condition there, contrary to what we did for the case $0 < \alpha < 1$. Nevertheless, we still impose $u|_{\Gamma_b} = 0$ in the numerical code, since the finite element formulation yields additional numerical dissipation that sufficiently regularizes the system.

In order to derive this $H^{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
therefore $\Psi_v \in W_{0,b}^{1,\gamma}(\text{BL})'$, which leads to

$$ \lim_{n \to \infty} \langle \Psi_v, u_{\epsilon_n} \rangle = \langle \Psi_v, u \rangle = \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z^\alpha \partial_z u \partial_z v. $$

Moreover,

$$ \left| \epsilon_n \int_{\text{BL}} \partial_z u_{\epsilon_n} \partial_z v \right| \le \epsilon_n \|u_{\epsilon_n}\|_{0;1,\gamma} \|v\|_{0;1,\gamma'} \le C(\gamma, G, f, \nu_h, \mu) \epsilon_n \|v\|_{0;1,\gamma'}, $$

giving

$$ \lim_{n \to \infty} \epsilon_n \int_{\text{BL}} \partial_z u_{\epsilon_n} \partial_z v = 0. $$

In addition, considering the properties of the sequence $(u_{\epsilon_n})_{n \in \mathbb{N}}$, we have

$$ \lim_{n \to \infty} \lambda \int_{\text{BL}} u_{\epsilon_n} v = \lambda \int_{\text{BL}} uv, \quad \lim_{n \to \infty} C_D \int_{\Gamma_{\text{top}}} u_{\epsilon_n} v = C_D \int_{\Gamma_{\text{top}}} uv. $$

Therefore $u$ satisfies (3.11).

It remains to show that the energy inequality (3.13) holds. Starting from (3.17), we have on one hand

$$ \lim_{n \to \infty} I_{1, \epsilon_n} = \lambda \int_{\text{BL}} |u|^2, \quad \lim_{n \to \infty} I_{3, \epsilon_n} = C_D \int_{\Gamma_{\text{top}}} |u|^2, $$

$$ \lim_{n \to \infty} I_{4, \epsilon_n} = \int_{\Gamma_{\text{top}}} u G, \quad \lim_{n \to \infty} I_{5, \epsilon_n} = \int_{\text{BL}} f u, $$

and on the other hand by lower semi-continuity

$$ \nu_h \int_{\text{BL}} |\nabla_h u|^2 + \mu \int_{\text{BL}} z^\alpha |\partial_z u|^2 \le \liminf_{n \to \infty} I_{2, \epsilon_n}, $$

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 $V_\alpha$, 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.

### 4. THE CASE $\alpha = 1$

We now consider the system:

$$ \begin{cases} \lambda u - \nu_h \Delta_h u - \mu \partial_z (z \partial_z u) = f & \text{in BL,} \\ C_D u + \mu z_{\text{top}} \frac{\partial u}{\partial z} = G & \text{on } \Gamma_{\text{top}}. \end{cases} \qquad (4.1) $$

As we shall see, the best we can get is an estimate in an $H^{1/2}$-like space near the bottom $\Gamma_b$. Therefore, we are not able to define the trace of such functions at $\Gamma_b$. This is why we do not impose any boundary condition there, contrary to what we did for the case $0 < \alpha < 1$. Nevertheless, we still impose $u|_{\Gamma_b} = 0$ in the numerical code, since the finite element formulation yields additional numerical dissipation that sufficiently regularizes the system.

In order to derive this $H^{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 $\mathcal{D}(\text{BL})$ denote the space of functions $u \in C^\infty(\text{BL})$ such that there exists $\delta = \delta(u) \in ]0, 1[$, ensuring that the support of $u$ is contained within $\Pi_2 \times [\delta(u), 1 - \delta(u)]$. This space is endowed with the standard Schwartz topology, characterized by the semi-norms $\|\hat{D}^n u\|_\infty$. In simpler terms, $\mathcal{D}(\text{BL})$ consists of $C^\infty$ functions that vanish in the vicinity of both the bottom and the top of the boundary layer. We use $\mathcal{D}'(\text{BL})$ to denote the corresponding space of distributions. Additionally, so far no risk of confusion arises, we still set

$$H_0^1(\text{BL}) = \{u \in H^1(\text{BL}), u = 0 \text{ on } \Gamma_b \cup \Gamma_{\text{top}}\}, \quad (4.2)$$

and $H^{-1}(\text{BL}) \subset \mathcal{D}'(\text{BL})$ is the topological dual of $H_0^1(\text{BL})$.

**Lemma 4.1.** Let $v \in \mathcal{D}'(\text{BL})$ such that $v \in H^{-1}(\text{BL})$ and also $\nabla v \in H^{-1}(\text{BL})$. Then $v \in L^2(\text{BL})$ and one has:

$$\|v\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})}. \quad (4.3)$$

The main consequence of Lemma 4.1 is the following.

**Corollary 4.1.** Let $V_1$ denotes the space:

$$V_1 = \{u \in \mathcal{D}'_b(\text{BL}) \text{ s.t. } u \in L^2(\text{BL}), \nabla_h u \in L^2(\text{BL}), z^{1/2} \partial_z u \in L^2(\text{BL})\}. \quad (4.4)$$

Then

$$V_1 \subset H^{1/2}(\text{BL}). \quad (4.5)$$

The space $V_1$ is the space that naturally suits problem (4.1) according to the energy balance (3.2) that still holds when $\alpha = 1$. It must be stressed that:

- showing that any solution of problem (4.1) is well-defined in $V_1$ by using (3.2) requires $\lambda > 0$;

- we have

$$V_1 \subset H^1(\Pi_2 \times [\delta, z_{\text{top}}]) \quad (4.6)$$

for any $\delta > 0$, which makes consistent the boundary condition at $\Gamma_{\text{top}}$.

**Definition 4.1.** We say that $u \in V_1$ is a very weak solution to (4.1) if for all $v \in V_1$,

$$\lambda \int_{\text{BL}} uv + \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z \partial_z u \partial_z v + \int_{\Gamma_{\text{top}}} (C_D u - G)v = \int_{\text{BL}} fv. \quad (4.7)$$

The main result of this section is the following.

**Theorem 4.1.** Let $f \in L^2$ and $G \in L^2$. Then, there exists a very weak solution $u \in V_1 \subset H^{1/2}(\text{BL})$ of problem (4.1), which satisfies

$$\|v\|_{V_1} \le C (\|f\|_{2,\text{BL}} + \|G\|_{2,\Gamma_{\text{top}}}).$$

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_h$-axis but not in the $z$-axis. This will be the purpose of the next sections.
utilize, let $\mathcal{D}(\text{BL})$ denote the space of functions $u \in C^\infty(\text{BL})$ such that there exists $\delta = \delta(u) \in ]0, 1[$, ensuring that the support of $u$ is contained within $\Pi_2 \times [\delta(u), 1 - \delta(u)]$. This space is endowed with the standard Schwartz topology, characterized by the semi-norms $\|\hat{D}^n u\|_\infty$. In simpler terms, $\mathcal{D}(\text{BL})$ consists of $C^\infty$ functions that vanish in the vicinity of both the bottom and the top of the boundary layer. We use $\mathcal{D}'(\text{BL})$ to denote the corresponding space of distributions. Additionally, so far no risk of confusion arises, we still set

$$H_0^1(\text{BL}) = \{u \in H^1(\text{BL}), u = 0 \text{ on } \Gamma_b \cup \Gamma_{\text{top}}\}, \quad (4.2)$$

and $H^{-1}(\text{BL}) \subset \mathcal{D}'(\text{BL})$ is the topological dual of $H_0^1(\text{BL})$.

**Lemma 4.1.** Let $v \in \mathcal{D}'(\text{BL})$ such that $v \in H^{-1}(\text{BL})$ and also $\nabla v \in H^{-1}(\text{BL})$. Then $v \in L^2(\text{BL})$ and one has:

$$\|v\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})}. \quad (4.3)$$

The main consequence of Lemma 4.1 is the following.

**Corollary 4.1.** Let $V_1$ denotes the space:

$$V_1 = \{u \in \mathcal{D}'_b(\text{BL}) \text{ s.t. } u \in L^2(\text{BL}), \nabla_h u \in L^2(\text{BL}), z^{1/2} \partial_z u \in L^2(\text{BL})\}. \quad (4.4)$$

Then

$$V_1 \subset H^{1/2}(\text{BL}). \quad (4.5)$$

The space $V_1$ is the space that naturally suits problem (4.1) according to the energy balance (3.2) that still holds when $\alpha = 1$. It must be stressed that:

- showing that any solution of problem (4.1) is well-defined in $V_1$ by using (3.2) requires $\lambda > 0$;

- we have

$$V_1 \subset H^1(\Pi_2 \times [\delta, z_{\text{top}}]) \quad (4.6)$$

for any $\delta > 0$, which makes consistent the boundary condition at $\Gamma_{\text{top}}$.

**Definition 4.1.** We say that $u \in V_1$ is a very weak solution to (4.1) if for all $v \in V_1$,

$$\lambda \int_{\text{BL}} uv + \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z \partial_z u \partial_z v + \int_{\Gamma_{\text{top}}} (C_D u - G)v = \int_{\text{BL}} fv. \quad (4.7)$$

The main result of this section is the following.

**Theorem 4.1.** Let $f \in L^2$ and $G \in L^2$. Then, there exists a very weak solution $u \in V_1 \subset H^{1/2}(\text{BL})$ of problem (4.1), which satisfies

$$\|v\|_{V_1} \le C (\|f\|_{2,\text{BL}} + \|G\|_{2,\Gamma_{\text{top}}}).$$

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_h$-axis but not in the $z$-axis. This will be the purpose of the next sections.
## 4.1. Framework

Let $\mathcal{T}_2 = \frac{2\pi}{L_x}\mathbb{Z} \times \frac{2\pi}{L_y}\mathbb{Z}$ and let $\mathbf{k} = (k_x, k_y) \in \mathcal{T}_2$ any wave vector. In the following, we set

$$ \mathcal{T}_{2,n} := \{ \mathbf{k} = (k_x, k_y) \in \mathcal{T}_2, |\mathbf{k}| \le 2\pi n \sqrt{1/L_x^2 + 1/L_y^2} \}. $$

Let $u \in L^1(\text{BL})$ and let $c_{\mathbf{k}} = c_{\mathbf{k}}(u; z)$ denotes the “horizontal” Fourier's coefficient at the wave vector $\mathbf{k}$, namely

$$ c_{\mathbf{k}}(u; z) := \frac{1}{\rho} \int_{\Pi_2} u(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h, $$

where $\rho = \sqrt{L_x^2 + L_y^2}$. Let $u_n$ be its partial sum of the Fourier series defined by

$$ u_n(\mathbf{x}_h, z) := \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}. \qquad (4.8) $$

A natural question is what assumptions are needed about $u$ to prove the convergence of the sequence $(u_n)_{n \in \mathbb{N}}$ toward $u$ (in some given topology). Before tackling this question, we consider the following spaces $\mathcal{D}(]0, z_{\text{top}}[) := C_0^{\infty}(]0, z_{\text{top}}[)$ and

$$ \mathcal{D}_n(\text{BL}) = \left\{ u = u_n(\mathbf{x}_h, z) = \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}, \ c_{\mathbf{k}} \in \mathcal{D}(]0, z_{\text{top}}[) \right\}, \qquad (4.9) $$

$$ V_n = \left\{ u = u_n(\mathbf{x}_h, z) = \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}, \ c_{\mathbf{k}} \in H_0^1(]0, z_{\text{top}}[) \right\}, \qquad (4.10) $$

$$ V_n^{-1} = \left\{ v = \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} \psi_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}, \ \psi_{\mathbf{k}} \in H^{-1}(]0, z_{\text{top}}[) \right\}. \qquad (4.11) $$

We will prove the following result.

**Lemma 4.2.** Let $u \in H_0^1(\text{BL})$ and $u_n$ be given by (4.8). Then, for all $n$, $u_n \in V_n$ and $(u_n)_{n \in \mathbb{N}}$ converges to $u$ in $H_0^1(\text{BL})$.

We first notice that the following isometries hold:

$$ \mathcal{D}_n(\text{BL}) \simeq \mathcal{D}(]0, z_{\text{top}}[)^n, \quad V_n \simeq H_0^1(]0, z_{\text{top}}[)^n, \quad V_n^{-1} \simeq H^{-1}([0, z_{\text{top}}])^n. \qquad (4.12) $$

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 $n \in \mathbb{N}$.

**Lemma 4.3.** Let $u_n \in \mathcal{D}'_n(\text{BL})$ such that $u_n, \nabla u_n \in V_n^{-1}$. Then $u_n \in L^2(\text{BL})$ and one has

$$ \|u_n\|_{0,2} \le \|u_n\|_{1,-1} + \|\nabla u_n\|_{1,-1}. \qquad (4.13) $$

It remains to pass to the limit in (4.13) when $n \to \infty$ to prove Lemma 4.1. We first must prove Lemma 4.2, which will be done step by step in the following subsections.
## 4.1. Framework

Let $\mathcal{T}_2 = \frac{2\pi}{L_x}\mathbb{Z} \times \frac{2\pi}{L_y}\mathbb{Z}$ and let $\mathbf{k} = (k_x, k_y) \in \mathcal{T}_2$ any wave vector. In the following, we set

$$ \mathcal{T}_{2,n} := \{ \mathbf{k} = (k_x, k_y) \in \mathcal{T}_2, |\mathbf{k}| \le 2\pi n \sqrt{1/L_x^2 + 1/L_y^2} \}. $$

Let $u \in L^1(\text{BL})$ and let $c_{\mathbf{k}} = c_{\mathbf{k}}(u; z)$ denotes the “horizontal” Fourier's coefficient at the wave vector $\mathbf{k}$, namely

$$ c_{\mathbf{k}}(u; z) := \frac{1}{\rho} \int_{\Pi_2} u(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h, $$

where $\rho = \sqrt{L_x^2 + L_y^2}$. Let $u_n$ be its partial sum of the Fourier series defined by

$$ u_n(\mathbf{x}_h, z) := \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}. \qquad (4.8) $$

A natural question is what assumptions are needed about $u$ to prove the convergence of the sequence $(u_n)_{n \in \mathbb{N}}$ toward $u$ (in some given topology). Before tackling this question, we consider the following spaces $\mathcal{D}(]0, z_{\text{top}}[) := C_0^{\infty}(]0, z_{\text{top}}[)$ and

$$ \mathcal{D}_n(\text{BL}) = \left\{ u = u_n(\mathbf{x}_h, z) = \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}, \ c_{\mathbf{k}} \in \mathcal{D}(]0, z_{\text{top}}[) \right\}, \qquad (4.9) $$

$$ V_n = \left\{ u = u_n(\mathbf{x}_h, z) = \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}, \ c_{\mathbf{k}} \in H_0^1(]0, z_{\text{top}}[) \right\}, \qquad (4.10) $$

$$ V_n^{-1} = \left\{ v = \sum_{\mathbf{k} \in \mathcal{T}_{2,n}} \psi_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}, \ \psi_{\mathbf{k}} \in H^{-1}(]0, z_{\text{top}}[) \right\}. \qquad (4.11) $$

We will prove the following result.

**Lemma 4.2.** Let $u \in H_0^1(\text{BL})$ and $u_n$ be given by (4.8). Then, for all $n$, $u_n \in V_n$ and $(u_n)_{n \in \mathbb{N}}$ converges to $u$ in $H_0^1(\text{BL})$.

We first notice that the following isometries hold:

$$ \mathcal{D}_n(\text{BL}) \simeq \mathcal{D}(]0, z_{\text{top}}[)^n, \quad V_n \simeq H_0^1(]0, z_{\text{top}}[)^n, \quad V_n^{-1} \simeq H^{-1}([0, z_{\text{top}}])^n. \qquad (4.12) $$

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 $n \in \mathbb{N}$.

**Lemma 4.3.** Let $u_n \in \mathcal{D}'_n(\text{BL})$ such that $u_n, \nabla u_n \in V_n^{-1}$. Then $u_n \in L^2(\text{BL})$ and one has

$$ \|u_n\|_{0,2} \le \|u_n\|_{1,-1} + \|\nabla u_n\|_{1,-1}. \qquad (4.13) $$

It remains to pass to the limit in (4.13) when $n \to \infty$ to prove Lemma 4.1. We first must prove Lemma 4.2, which will be done step by step in the following subsections.
## 4.2. $L^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 $u \in H^2(\text{BL})$, $u_n$ its Fourier's series expansion as given by (4.8). Then $u_n \to u$ in $L^2(\text{BL})$ as $n \to \infty$.

**Proof.** Step 1: Estimate of the $c_k$. Let $k \in T_{2,n} \setminus \{0\}$. By Fubini's Theorem and two integration by parts, we get

$$ 
\begin{aligned}
c_k(u; z) &= \frac{1}{\rho k_x^2} \int_{\Pi_2} \frac{\partial^2 u}{\partial x^2}(\mathbf{x}_h, z) e^{-ik \cdot \mathbf{x}_h} \mathrm{d}\mathbf{x}_h, && \text{if } k_x \neq 0, \\
c_k(u; z) &= -\frac{1}{\rho k_y^2} \int_{\Pi_2} \frac{\partial^2 u}{\partial y^2}(\mathbf{x}_h, z) e^{-ik \cdot \mathbf{x}_h} \mathrm{d}\mathbf{x}_h, && \text{if } k_y \neq 0.
\end{aligned} 
\quad (4.14) $$

Consequently, we get

$$ |c_k(u; z)| \le \frac{1}{\rho} \inf \left\{ \frac{1}{k_x^2}, \frac{1}{k_y^2} \right\} \int_{\Pi_2} |\nabla_h u(\mathbf{x}_h, z)|^2 \mathrm{d}\mathbf{x}_h \le \frac{C}{\rho |\mathbf{k}|^2} \|\nabla_h u(\cdot, z)\|_{L^2(\Pi_2)}^2. \quad (4.15) $$

Therefore, as $u \in H^1(\text{BL})$, $z \to c_k(u; z) \in L^2([0, z_{\text{top}}])$ and is finite almost everywhere.

**Step 2:** Convergence. We deduce from classical results that for almost all $z \in ]0, z_{\text{top}}[$,

$$ u_n(\cdot, z) \to u(\cdot, z) \quad \text{in } L^2(\Pi_2). \quad (4.16) $$

Put another way:

$$ \varepsilon_n(z) = \int_{\Pi_2} |u_n(\mathbf{x}_h, z) - u(\mathbf{x}_h, z)|^2 \mathrm{d}\mathbf{x}_h \to 0 \quad \text{as } n \to \infty. \quad (4.17) $$

Moreover, by the previous step,

$$ 0 \le \varepsilon_n \le \sum_{|\mathbf{k}| \ge \lambda n} |c_k(u; z)|^2 \le \frac{C}{\rho} \|\nabla_h u\|_{L^2(\Pi_2)}^2 \left( \sum_{|\mathbf{k}| \ge \lambda n} \frac{1}{|\mathbf{k}|^4} \right) = R_n \|\nabla_h u\|_{L^2(\Pi_2)}^2, \quad (4.18) $$

giving

$$ \|u - u_n\|_{0,2}^2 = \int_0^{z_{\text{top}}} \varepsilon_n(z) \mathrm{d}z \le R_n \|\nabla_h u\|_{0,2}^2 \to 0 \quad \text{as } n \to \infty, \quad (4.19) $$

which concludes the proof.

## 4.3. Differentiability of the coefficient $c_k(u; z)$

**Lemma 4.5.** Let $u \in H^1(\text{BL})$ such that $\frac{\partial u}{\partial z} \in H^1(\text{BL})$. Then $c_k(u; \cdot) \in H^1([0, z_{\text{top}}])$ and one has for almost all $z \in ]0, z_{\text{top}}[$,

$$ \frac{\mathrm{d}}{\mathrm{d}z} c_k(u; z) = \int_{\Pi_2} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) e^{-ik \cdot \mathbf{x}_h} \mathrm{d}\mathbf{x}_h. \quad (4.20) $$

**Proof.** Let us write

$$ \frac{\partial u}{\partial z}(\mathbf{x}_h, z) = \frac{\partial u}{\partial z}(\mathbf{x}_h, 0) + \int_0^{z_{\text{top}}} \frac{\partial^2 u}{\partial z^2}(\mathbf{x}_h, z') \mathrm{d}z'. \quad (4.21) $$
## 4.2. $L^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 $u \in H^2(\text{BL})$, $u_n$ its Fourier's series expansion as given by (4.8). Then $u_n \to u$ in $L^2(\text{BL})$ as $n \to \infty$.

**Proof.** Step 1: Estimate of the $c_k$. Let $k \in T_{2,n} \setminus \{0\}$. By Fubini's Theorem and two integration by parts, we get

$$ 
\begin{aligned}
c_k(u; z) &= \frac{1}{\rho k_x^2} \int_{\Pi_2} \frac{\partial^2 u}{\partial x^2}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} \mathrm{d}\mathbf{x}_h, && \text{if } k_x \neq 0, \\
c_k(u; z) &= -\frac{1}{\rho k_y^2} \int_{\Pi_2} \frac{\partial^2 u}{\partial y^2}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} \mathrm{d}\mathbf{x}_h, && \text{if } k_y \neq 0.
\end{aligned} 
\quad (4.14) $$

Consequently, we get

$$ |c_k(u; z)| \le \frac{1}{\rho} \inf \left\{ \frac{1}{k_x^2}, \frac{1}{k_y^2} \right\} \int_{\Pi_2} |\nabla_h u(\mathbf{x}_h, z)|^2 \mathrm{d}\mathbf{x}_h \le \frac{C}{\rho |\mathbf{k}|^2} \|\nabla_h u(\cdot, z)\|_{L^2(\Pi_2)}^2. \quad (4.15) $$

Therefore, as $u \in H^1(\text{BL})$, $z \to c_k(u; z) \in L^2([0, z_{\text{top}}])$ and is finite almost everywhere.

**Step 2:** Convergence. We deduce from classical results that for almost all $z \in ]0, z_{\text{top}}[$,

$$ u_n(\cdot, z) \to u(\cdot, z) \quad \text{in } L^2(\Pi_2). \quad (4.16) $$

Put another way:

$$ \varepsilon_n(z) = \int_{\Pi_2} |u_n(\mathbf{x}_h, z) - u(\mathbf{x}_h, z)|^2 \mathrm{d}\mathbf{x}_h \to 0 \quad \text{as } n \to \infty. \quad (4.17) $$

Moreover, by the previous step,

$$ 0 \le \varepsilon_n \le \sum_{|\mathbf{k}| \ge \lambda n} |c_k(u; z)|^2 \le \frac{C}{\rho} \|\nabla_h u\|_{L^2(\Pi_2)}^2 \left( \sum_{|\mathbf{k}| \ge \lambda n} \frac{1}{|\mathbf{k}|^4} \right) = R_n \|\nabla_h u\|_{L^2(\Pi_2)}^2, \quad (4.18) $$

giving

$$ \|u - u_n\|_{0,2}^2 = \int_0^{z_{\text{top}}} \varepsilon_n(z) \mathrm{d}z \le R_n \|\nabla_h u\|_{0,2}^2 \to 0 \quad \text{as } n \to \infty, \quad (4.19) $$

which concludes the proof.

## 4.3. Differentiability of the coefficient $c_k(u; z)$

**Lemma 4.5.** Let $u \in H^1(\text{BL})$ such that $\frac{\partial u}{\partial z} \in H^1(\text{BL})$. Then $c_k(u; \cdot) \in H^1([0, z_{\text{top}}])$ and one has for almost all $z \in ]0, z_{\text{top}}[$,

$$ \frac{\mathrm{d}}{\mathrm{d}z} c_k(u; z) = \int_{\Pi_2} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} \mathrm{d}\mathbf{x}_h. \quad (4.20) $$

**Proof.** Let us write

$$ \frac{\partial u}{\partial z}(\mathbf{x}_h, z) = \frac{\partial u}{\partial z}(\mathbf{x}_h, 0) + \int_0^{z_{\text{top}}} \frac{\partial^2 u}{\partial z^2}(\mathbf{x}_h, z') \mathrm{d}z'. \quad (4.21) $$
By the trace theorem,

$$ \frac{\partial u}{\partial z}(\mathbf{x}_h, 0) \in H^{\frac{1}{2}}(\Gamma_b) \hookrightarrow L^2(\Gamma_b) \simeq L^2(\Pi_2) $$

hence, $\forall z \in ]0, z_{\text{top}}[$,

$$ \left| \frac{\partial u}{\partial z}(\mathbf{x}_h, z) \right| \le \left| \frac{\partial u}{\partial z}(\mathbf{x}_h, 0) \right| + \int_0^{z_{\text{top}}} \left| \frac{\partial^2 u}{\partial z^2}(\mathbf{x}_h, z') \right| dz' \in L^1(\Pi_2). $$

Therefore, formula (4.20) is a classical consequence of the Lebesgue Theorem. □

**Lemma 4.6.** Let $u \in H_0^1(\text{BL})$. Then $c_k(u; \cdot) \in H_0^1(]0, z_{\text{top}}[)$ and (4.20) holds a.e in $]0, z_{\text{top}}[$.

**Proof.** Step 1: Approximations. Let $\varepsilon > 0$, and $u_\varepsilon \in H_0^1(\text{BL}) \cap H^2(\text{BL})$ such that

$$ -\varepsilon \Delta u_\varepsilon + u_\varepsilon = u. \qquad (4.22) $$

According to standard results about the Helmholtz equation, we know that that $u_\varepsilon \to u$ as $\varepsilon \to 0$, strongly in $L^2(\text{BL})$, weakly in $H_0^1(\text{BL})$. Moreover, by Lemma 4.5, we also know that

$$ \frac{d}{dz} c_k(u_\varepsilon; z) = \int_{\Pi_2} \frac{\partial u_\varepsilon}{\partial z}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h \qquad (4.23) $$

as well as $c_k(u_\varepsilon; z) \to c_k(u; z)$ in $L^2(]0, z_{\text{top}}[)$.

**Step 2:** Derivative in the sense of the distributions. Let $\varphi \in \mathcal{D}(]0, z_{\text{top}}[)$. The starting point is the identity

$$ - \int_{0}^{z_{\text{top}}} \frac{d\varphi}{dz}(z) c_{\mathbf{k}}(u_{\varepsilon}; z) dz = \int_{0}^{z_{\text{top}}} \frac{d}{dz} c_{\mathbf{k}}(u_{\varepsilon}; z) \varphi(z) dz. \qquad (4.24) $$

From the results of step 1, we already know that

$$ \int_{0}^{z_{\text{top}}} \frac{d\varphi}{dz}(z) c_{\mathbf{k}}(u_{\varepsilon}; z) dz \xrightarrow{\varepsilon \to 0} \int_{0}^{z_{\text{top}}} \frac{d\varphi}{dz}(z) c_{\mathbf{k}}(u; z) dz. \qquad (4.25) $$

We must pass to the limit in the r.h.s of (4.24). By (4.20), we have

$$ \int_{0}^{z_{\text{top}}} \frac{d}{dz} c_{\mathbf{k}}(u_{\varepsilon}; z) \varphi(z) dz = \int_{\text{BL}} \frac{\partial u_{\varepsilon}}{\partial z}(\mathbf{x}_h, z) \varphi(z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz. $$

Therefore, as $u_\varepsilon \to u$ weakly in $H_0^1(\text{BL})$,

$$ \int_{\text{BL}} \frac{\partial u_\varepsilon}{\partial z}(\mathbf{x}_h, z) \varphi(z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz \xrightarrow{\varepsilon \to 0} \int_{\text{BL}} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) \varphi(z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz \\ = \int_{0}^{z_{\text{top}}} \varphi(z) \int_{\Pi_2} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz. \qquad (4.26) $$

By combining (4.24)–(4.26), we see that (4.20) holds in the sense of the distribution, which concludes the proof because

$$ z \to \int_{\Pi_2} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) e^{i-\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h \in L^2(]0, z_{\text{top}}[). $$

By using the same argument, we also have the following, where we note $P_n u = u_n$.

**Lemma 4.7.** Let $u \in L^2(\text{BL})$. Then $P_n u$ converges to $u$ in $L^2(\text{BL})$ as $n \to \infty$.
By the trace theorem,

$$ \frac{\partial u}{\partial z}(\mathbf{x}_h, 0) \in H^{\frac{1}{2}}(\Gamma_b) \hookrightarrow L^2(\Gamma_b) \simeq L^2(\Pi_2) $$

hence, $\forall z \in ]0, z_{\text{top}}[$,

$$ \left| \frac{\partial u}{\partial z}(\mathbf{x}_h, z) \right| \le \left| \frac{\partial u}{\partial z}(\mathbf{x}_h, 0) \right| + \int_0^{z_{\text{top}}} \left| \frac{\partial^2 u}{\partial z^2}(\mathbf{x}_h, z') \right| dz' \in L^1(\Pi_2). $$

Therefore, formula (4.20) is a classical consequence of the Lebesgue Theorem. □

**Lemma 4.6.** Let $u \in H_0^1(\text{BL})$. Then $c_k(u; \cdot) \in H_0^1(]0, z_{\text{top}}[)$ and (4.20) holds a.e in $]0, z_{\text{top}}[$.

**Proof.** Step 1: Approximations. Let $\varepsilon > 0$, and $u_\varepsilon \in H_0^1(\text{BL}) \cap H^2(\text{BL})$ such that

$$ -\varepsilon \Delta u_\varepsilon + u_\varepsilon = u. \qquad (4.22) $$

According to standard results about the Helmholtz equation, we know that that $u_\varepsilon \to u$ as $\varepsilon \to 0$, strongly in $L^2(\text{BL})$, weakly in $H_0^1(\text{BL})$. Moreover, by Lemma 4.5, we also know that

$$ \frac{d}{dz} c_k(u_\varepsilon; z) = \int_{\Pi_2} \frac{\partial u_\varepsilon}{\partial z}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h \qquad (4.23) $$

as well as $c_k(u_\varepsilon; z) \to c_k(u; z)$ in $L^2(]0, z_{\text{top}}[)$.

**Step 2:** Derivative in the sense of the distributions. Let $\varphi \in \mathcal{D}(]0, z_{\text{top}}[)$. The starting point is the identity

$$ -\int_{0}^{z_{\text{top}}} \frac{d\varphi}{dz}(z) c_k(u_{\varepsilon}; z) dz = \int_{0}^{z_{\text{top}}} \frac{d}{dz} c_k(u_{\varepsilon}; z) \varphi(z) dz. \qquad (4.24) $$

From the results of step 1, we already know that

$$ \int_{0}^{z_{\text{top}}} \frac{d\varphi}{dz}(z) c_k(u_{\varepsilon}; z) dz \xrightarrow{\varepsilon \to 0} \int_{0}^{z_{\text{top}}} \frac{d\varphi}{dz}(z) c_k(u; z) dz. \qquad (4.25) $$

We must pass to the limit in the r.h.s of (4.24). By (4.20), we have

$$ \int_{0}^{z_{\text{top}}} \frac{d}{dz} c_k(u_{\varepsilon}; z) \varphi(z) dz = \int_{\text{BL}} \frac{\partial u_{\varepsilon}}{\partial z}(\mathbf{x}_h, z) \varphi(z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz. $$

Therefore, as $u_\varepsilon \to u$ weakly in $H_0^1(\text{BL})$,

$$ \int_{\text{BL}} \frac{\partial u_{\varepsilon}}{\partial z}(\mathbf{x}_h, z) \varphi(z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz \xrightarrow{\varepsilon \to 0} \int_{\text{BL}} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) \varphi(z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz \\ = \int_{0}^{z_{\text{top}}} \varphi(z) \int_{\Pi_2} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) e^{-i\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h dz. \qquad (4.26) $$

By combining (4.24)–(4.26), we see that (4.20) holds in the sense of the distribution, which concludes the proof because

$$ z \to \int_{\Pi_2} \frac{\partial u}{\partial z}(\mathbf{x}_h, z) e^{i-\mathbf{k} \cdot \mathbf{x}_h} d\mathbf{x}_h \in L^2(]0, z_{\text{top}}[). $$

By using the same argument, we also have the following, where we note $P_n u = u_n$.

**Lemma 4.7.** Let $u \in L^2(\text{BL})$. Then $P_n u$ converges to $u$ in $L^2(\text{BL})$ as $n \to \infty$.
Proof. Let $u \in L^2(\text{BL})$ and $u_\varepsilon \in H^2(\text{BL})$ given by (4.22). Then, we know from Lemma 4.4 that for any $\varepsilon > 0$ be fixed, $P_n u_\varepsilon$ converges to $u_\varepsilon$ in $L^2(\text{BL})$. We write

$$P_n u - u = P_n(u - u_\varepsilon) + (P_n u_\varepsilon - u_\varepsilon) + (u - u_\varepsilon).$$

Then, since $\|P_n\| \le 1$, we get

$$\|P_n u - u\|_{0,2} \le 2\|u - u_\varepsilon\|_{0,2} + \|P_n u_\varepsilon - u_\varepsilon\|_{0,2}.$$

The rest of the proof is straightforward. □

## 4.4. $H^1$ Convergence

We are now able to prove Lemma 4.2, stating the convergence of the sequence $(u_n)_{n \in \mathbb{N}}$ to $u$ in $H^1(\text{BL})$. We already know from Lemma 4.7 that $u_n = P_n u \underset{n \to \infty}{\longrightarrow} u$ in $L^2(\text{BL})$. Let $\mathcal{F}_n = \mathcal{F}_n(\mathbf{x}_h)$ denotes the Fejer kernel over $\Pi_2$. Then, for a.e $z \in ]0, z_{\text{top}}[$, $u_n$ is given by

$$u_n(\mathbf{x}_h, z) = \int_{\Pi_2} \mathcal{F}_n(\mathbf{y}_h)u(\mathbf{x}_h + \mathbf{y}_h, z) \, d\mathbf{y}_h = \mathcal{F}_n(\cdot) \star u(\cdot, z)(\mathbf{x}_h). \quad (4.27)$$

As the Lebesgue measure over BL is $\sigma$-finite, we deduce from the Lebesgue–Fubini theorem that for a.e $z \in ]0, z_{\text{top}}[$, $u(\cdot, z) \in H^1(\Pi_2)$. Therefore, always a.e $z \in ]0, z_{\text{top}}[$,

$$\nabla_h u_n(\mathbf{x}_h, z) = \mathcal{F}_n(\cdot) \star \nabla_h u(\cdot, z)(\mathbf{x}_h), \quad (4.28)$$

and by standard results,

$$\varepsilon_n(z) = \int_{\Pi_2} |\nabla_h u_n(\mathbf{x}_h, z) - \nabla_h u(\mathbf{x}_h, z)|^2 d\mathbf{x}_h \underset{n \to \infty}{\longrightarrow} 0. \quad (4.29)$$

Moreover, by the Young inequality,

$$\begin{aligned} 0 \le \varepsilon_n(z) &\le (\|\mathcal{F}_n\|_{\Pi_2;0,1} \|\nabla_h u(\cdot, z)\|_{\Pi_2;0,2} + \|\nabla_h u(\cdot, z)\|_{\Pi_2;0,2})^2 \\ &\le 4 \|\nabla_h u(\cdot, z)\|_{\Pi_2;0,2}^2 \in L^1(]0, z_{\text{top}}[), \end{aligned} \quad (4.30)$$

where $\|\cdot\|_{\Pi_2; s, p}$ denotes the usual norm of the Sobolev space $W^{s,p}(\Pi_2)$. Therefore, by the Lebesgue Theorem, we get

$$\int_{0}^{z_{\text{top}}} \varepsilon_n(z) = \int_{\text{BL}} |u_n(\mathbf{x}_h, z) - u(\mathbf{x}_h, z)|^2 d\mathbf{x}_h \, dz \underset{n \to \infty}{\longrightarrow} 0, \quad (4.31)$$

which, put in another way, yields $\nabla_h u_n \underset{n \to \infty}{\longrightarrow} \nabla_h u$ in $L^2(\text{BL})$. Similarly, we deduce from Lemma 4.6 that

$$\frac{\partial u_n}{\partial z}(\mathbf{x}_h, z) = \mathcal{F}_n(\cdot) \star \frac{\partial u}{\partial z}(\cdot, z)(\mathbf{x}_h), \quad (4.32)$$

and by an analysis similar to the previous one, we can conclude that

$$\frac{\partial u_n}{\partial z} \underset{n \to \infty}{\longrightarrow} \frac{\partial u}{\partial z},$$

in $L^2(\text{BL})$, which finishes the proof of Lemma 4.2.
Proof. Let $u \in L^2(\text{BL})$ and $u_\varepsilon \in H^2(\text{BL})$ given by (4.22). Then, we know from Lemma 4.4 that for any $\varepsilon > 0$ be fixed, $P_n u_\varepsilon$ converges to $u_\varepsilon$ in $L^2(\text{BL})$. We write

$$P_n u - u = P_n(u - u_\varepsilon) + (P_n u_\varepsilon - u_\varepsilon) + (u - u_\varepsilon).$$

Then, since $\|P_n\| \le 1$, we get

$$\|P_n u - u\|_{0,2} \le 2\|u - u_\varepsilon\|_{0,2} + \|P_n u_\varepsilon - u_\varepsilon\|_{0,2}.$$

The rest of the proof is straightforward. □

## 4.4. $H^1$ Convergence

We are now able to prove Lemma 4.2, stating the convergence of the sequence $(u_n)_{n \in \mathbb{N}}$ to $u$ in $H^1(\text{BL})$. We already know from Lemma 4.7 that $u_n = P_n u \underset{n \to \infty}{\longrightarrow} u$ in $L^2(\text{BL})$. Let $\mathcal{F}_n = \mathcal{F}_n(\mathbf{x}_h)$ denotes the Fejer kernel over $\Pi_2$. Then, for a.e $z \in ]0, z_{\text{top}}[$, $u_n$ is given by

$$u_n(\mathbf{x}_h, z) = \int_{\Pi_2} \mathcal{F}_n(\mathbf{y}_h)u(\mathbf{x}_h + \mathbf{y}_h, z) \, d\mathbf{y}_h = \mathcal{F}_n(\cdot) \star u(\cdot, z)(\mathbf{x}_h). \quad (4.27)$$

As the Lebesgue measure over BL is $\sigma$-finite, we deduce from the Lebesgue–Fubini theorem that for a.e $z \in ]0, z_{\text{top}}[$, $u(\cdot, z) \in H^1(\Pi_2)$. Therefore, always a.e $z \in ]0, z_{\text{top}}[$,

$$\nabla_h u_n(\mathbf{x}_h, z) = \mathcal{F}_n(\cdot) \star \nabla_h u(\cdot, z)(\mathbf{x}_h), \quad (4.28)$$

and by standard results,

$$\varepsilon_n(z) = \int_{\Pi_2} |\nabla_h u_n(\mathbf{x}_h, z) - \nabla_h u(\mathbf{x}_h, z)|^2 \, d\mathbf{x}_h \underset{n \to \infty}{\longrightarrow} 0. \quad (4.29)$$

Moreover, by the Young inequality,

$$\begin{aligned} 0 \le \varepsilon_n(z) &\le (\|\mathcal{F}_n\|_{\Pi_2;0,1} \|\nabla_h u(\cdot, z)\|_{\Pi_2;0,2} + \|\nabla_h u(\cdot, z)\|_{\Pi_2;0,2})^2 \\ &\le 4 \|\nabla_h u(\cdot, z)\|_{\Pi_2;0,2}^2 \in L^1(]0, z_{\text{top}}[), \end{aligned} \quad (4.30)$$

where $\|\cdot\|_{\Pi_2; s, p}$ denotes the usual norm of the Sobolev space $W^{s, p}(\Pi_2)$. Therefore, by the Lebesgue Theorem, we get

$$\int_{0}^{z_{\text{top}}} \varepsilon_n(z) = \int_{\text{BL}} |u_n(\mathbf{x}_h, z) - u(\mathbf{x}_h, z)|^2 \, d\mathbf{x}_h \, dz \underset{n \to \infty}{\longrightarrow} 0, \quad (4.31)$$

which, put in another way, yields $\nabla_h u_n \underset{n \to \infty}{\longrightarrow} \nabla_h u$ in $L^2(\text{BL})$. Similarly, we deduce from Lemma 4.6 that

$$\frac{\partial u_n}{\partial z}(\mathbf{x}_h, z) = \mathcal{F}_n(\cdot) \star \frac{\partial u}{\partial z}(\cdot, z)(\mathbf{x}_h), \quad (4.32)$$

and by an analysis similar to the previous one, we can conclude that

$$\frac{\partial u_n}{\partial z} \underset{n \to \infty}{\longrightarrow} \frac{\partial u}{\partial z},$$

in $L^2(\text{BL})$, which finishes the proof of Lemma 4.2.
### 4.5. Interpolation spaces

As a consequence of Lemmas 4.2 and 4.7, combined with the interpolation theorem proved by Lions and Magenes [12], we have the following.

**Lemma 4.8.** Let $s \in [0, 1]$, and let

$$H^s(\text{BL}) = [L^2(\text{BL}), H^1(\text{BL})]_s = \mathbb{D}(\nabla^s).$$

Then, $u = \sum_{\mathbf{k} \in \mathcal{T}_2} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h} \in H^s(\text{BL})$ if and only if $\forall \mathbf{k} \in \mathcal{T}_2$, $c_{\mathbf{k}} \in H^s([0, z_{\text{top}}])$ and

$$\|u\|_{s,2}^2 = \sum_{\mathbf{k} \in \mathcal{T}_2} \left( |\mathbf{k}|^{2s} \int_0^{z_{\text{top}}} |c_{\mathbf{k}}(z)|^2 \, dz + \|c_{\mathbf{k}}\|_{s,2}^2 \right) < \infty.$$

We now need the definition of the space $H_{00}^{1/2}(\text{BL})$ suited to our geometry, and its dual space. Let $\rho = \rho(z)$ be a $C^\infty$ non negative function on $[0, z_{\text{top}}]$, and such that

$$\lim_{z \to 0} \frac{\rho(z)}{z} = \lim_{z \to z_{\text{top}}} \frac{\rho(z)}{z_{\text{top}} - z} = 1.$$

Then in our case,

$$H_{00}^{1/2}(\text{BL}) = \{u \in H^{1/2}(\text{BL}) \text{ s.t. } \rho^{-1/2}u \in L^2(\text{BL})\},$$

equipped with the norm

$$\|u\|_{H_{00}^{1/2}} = \left( \|u\|_{1/2,2}^2 + \|\rho^{-1/2}u\|_{0,2}^2 \right)^{\frac{1}{2}}.$$

According to Lemma 4.8, we have the following.

**Lemma 4.9.** Let $u = \sum_{\mathbf{k} \in \mathcal{T}_2} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}$. Then, $u \in H_{00}^{1/2}(\text{BL})$ if and only if $\forall \mathbf{k} \in \mathcal{T}_2$, $c_{\mathbf{k}} \in H_{00}^{1/2}([0, z_{\text{top}}])$ and

$$\|u\|_{H_{00}^{1/2}}^2 = \|u\|_{1/2,2}^2 + \sum_{\mathbf{k} \in \mathcal{T}_2} \int_0^{z_{\text{top}}} \frac{|c_{\mathbf{k}}(z)|^2}{\rho(z)} \, dz < \infty.$$

Now following [12], let us consider the linear operator $\Lambda = \nabla^{-1} : L^2 \to H^{-1}$, defined such that

$$\forall u, v \in L^2, \quad (u, v)_2 = (\Lambda u, \Lambda v)_{-1}, \qquad (4.33)$$

with natural notations to denote the various scalar products. Then, the space $[H^{-1}, L^2]_{1/2}$ is the domain of $\Lambda^{1/2}$, and one has, according to the above,

$$[H^{-1}, L^2]_{1/2} = (H_{00}^{1/2}(\text{BL}))'. \qquad (4.34)$$

In order to complete the interpolation tool box, let us consider $W^{-1}$ as the closure of $C^\infty(\text{BL})$ subjected to the norm

$$\|v\|_{W^{-1}} = \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})},$$

$W^0$ is the closure with respect to the norm

$$\|\nabla v\|_{L^2(\text{BL})} + \|v\|_{L^2(\text{BL})}.$$

Note that $W^0 \approx H^1(\text{BL})$. The question is the characterization of the interpolation space $[W^{-1}, W^0]_{1/2}$. From (4.33), we have

$$\forall u, v \in W^0, \quad (u, v)_2 + (\nabla u, \nabla v)_2 = (\Lambda u, \Lambda v)_{-1} + (\Lambda \nabla u, \Lambda \nabla v)_{-1},$$

which allows to characterize the space $[W^{-1}, W^0]_{1/2}$ thanks to (4.34) by its norm given by:

$$\|v\|_{[W^{-1}, W^0]_{1/2}} = \|v\|_{(H_{00}^{1/2}(\text{BL}))'} + \|\nabla v\|_{(H_{00}^{1/2}(\text{BL}))'}. \qquad (4.35)$$
### 4.5. Interpolation spaces

As a consequence of Lemmas 4.2 and 4.7, combined with the interpolation theorem proved by Lions and Magenes [12], we have the following.

**Lemma 4.8.** Let $s \in [0, 1]$, and let

$$H^s(\text{BL}) = [L^2(\text{BL}), H^1(\text{BL})]_s = \mathbb{D}(\nabla^s).$$

Then, $u = \sum_{\mathbf{k} \in \mathcal{T}_2} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h} \in H^s(\text{BL})$ if and only if $\forall \mathbf{k} \in \mathcal{T}_2$, $c_{\mathbf{k}} \in H^s([0, z_{\text{top}}])$ and

$$\|u\|_{s,2}^2 = \sum_{\mathbf{k} \in \mathcal{T}_2} \left( |\mathbf{k}|^{2s} \int_0^{z_{\text{top}}} |c_{\mathbf{k}}(z)|^2 \, dz + \|c_{\mathbf{k}}\|_{s,2}^2 \right) < \infty.$$

We now need the definition of the space $H_{00}^{1/2}(\text{BL})$ suited to our geometry, and its dual space. Let $\rho = \rho(z)$ be a $C^\infty$ non negative function on $[0, z_{\text{top}}]$, and such that

$$\lim_{z \to 0} \frac{\rho(z)}{z} = \lim_{z \to z_{\text{top}}} \frac{\rho(z)}{z_{\text{top}} - z} = 1.$$

Then in our case,

$$H_{00}^{1/2}(\text{BL}) = \{u \in H^{1/2}(\text{BL}) \text{ s.t. } \rho^{-1/2}u \in L^2(\text{BL})\},$$

equipped with the norm

$$\|u\|_{H_{00}^{1/2}} = \left( \|u\|_{1/2,2}^2 + \|\rho^{-1/2}u\|_{0,2}^2 \right)^{\frac{1}{2}}.$$

According to Lemma 4.8, we have the following.

**Lemma 4.9.** Let $u = \sum_{\mathbf{k} \in \mathcal{T}_2} c_{\mathbf{k}}(z) e^{i\mathbf{k} \cdot \mathbf{x}_h}$. Then, $u \in H_{00}^{1/2}(\text{BL})$ if and only if $\forall \mathbf{k} \in \mathcal{T}_2$, $c_{\mathbf{k}} \in H_{00}^{1/2}([0, z_{\text{top}}])$ and

$$\|u\|_{H_{00}^{1/2}}^2 = \|u\|_{1/2,2}^2 + \sum_{\mathbf{k} \in \mathcal{T}_2} \int_0^{z_{\text{top}}} \frac{|c_{\mathbf{k}}(z)|^2}{\rho(z)} \, dz < \infty.$$

Now following [12], let us consider the linear operator $\Lambda = \nabla^{-1} : L^2 \to H^{-1}$, defined such that

$$\forall u, v \in L^2, \quad (u, v)_2 = (\Lambda u, \Lambda v)_{-1}, \qquad (4.33)$$

with natural notations to denote the various scalar products. Then, the space $[H^{-1}, L^2]_{1/2}$ is the domain of $\Lambda^{1/2}$, and one has, according to the above,

$$[H^{-1}, L^2]_{1/2} = (H_{00}^{1/2}(\text{BL}))'. \qquad (4.34)$$

In order to complete the interpolation tool box, let us consider $W^{-1}$ as the closure of $C^\infty(\text{BL})$ subjected to the norm

$$\|v\|_{W^{-1}} = \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})},$$

$W^0$ is the closure with respect to the norm

$$\|\nabla v\|_{L^2(\text{BL})} + \|v\|_{L^2(\text{BL})}.$$

Note that $W^0 \approx H^1(\text{BL})$. The question is the characterization of the interpolation space $[W^{-1}, W^0]_{1/2}$. From (4.33), we have

$$\forall u, v \in W^0, \quad (u, v)_2 + (\nabla u, \nabla v)_2 = (\Lambda u, \Lambda v)_{-1} + (\Lambda \nabla u, \Lambda \nabla v)_{-1},$$

which allows to characterize the space $[W^{-1}, W^0]_{1/2}$ thanks to (4.34) by its norm given by:

$$\|v\|_{[W^{-1}, W^0]_{1/2}} = \|v\|_{(H_{00}^{1/2}(\text{BL}))'} + \|\nabla v\|_{(H_{00}^{1/2}(\text{BL}))'}. \qquad (4.35)$$
## 4.6. 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,

$$ \|v\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})} \quad (4.36) $$

satisfied by any distribution $v \in D'(\text{BL})$ such that $v, \nabla v \in H^{-1}(\text{BL})$, the dual of the space $H_0^1(\text{BL})$ as defined by (4.2). Then, we will prove Corollary 4.1, namely

$$ V_1 \subset H^{\frac{1}{2}}(\text{BL}), \quad (4.37) $$

with continuous injection, where $V_1$ is defined by (4.4), equipped with the norm

$$ \|u\|_{V^1} = \|u\|_{L^2(\text{BL})} + \|\nabla_h u\|_{L^2(\text{BL})} + \|z^{1/2} \partial_z u\|_{L^2(\text{BL})}. $$

In particular, being $V_1$ a closed subspace of $H^{1/2}(\text{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 $v \in D'(\text{BL})$ such that $v, \nabla v \in H^{-1}$, and $v_n, \nabla v_n \in V_n^{-1}$ be given by

$$ \forall u \in V_n, U \in V_n^3, \quad \langle v_n, u \rangle = \langle v, u \rangle, \quad \langle \nabla v_n, U \rangle = \langle \nabla v, U \rangle, $$

where $V_n^{-1}$ is defined by (4.11). We have, for all $n \in \mathbb{N}$,

$$ \|\nabla v_n\|_{H^{-1}(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})}, \quad \|v_n\|_{H^{-1}(\text{BL})} \le \|v\|_{H^{-1}(\text{BL})}, $$

and note that (by explicit computation) it easily follows that the constants on the right-hand side are equal to 1, independently of $n \in \mathbb{N}$. Therefore, according to inequality (4.13),

$$ \|v_n\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})}. \quad (4.38) $$

Then, the sequence $(v_n)_{n \in \mathbb{N}}$ is bounded in $L^2(\text{BL})$, and we can extract a subsequence (still denoted by $(v_n)_{n \in \mathbb{N}}$) that weakly converges to some $\tilde{v}$, which is equal to $v$ in $H^{-1}(\text{BL})$, therefore in $L^2(\text{BL})$, which yields $v \in L^2(\text{BL})$ and

$$ \|v\|_{L^2(\text{BL})} \le \liminf_{n \to \infty} \|v_n\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})}, $$

hence (4.36).

*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

$$ \mathrm{Id} : W^{-1} \to L^2 $$

is continuous. Moreover, $\mathrm{Id} : W^0 \to H^1$ is also continuous. Therefore, by the interpolation theorem [12], also the restriction of the identity (denoted still by $\mathrm{Id}$)

$$ \mathrm{Id} : [W^{-1}, W^0]_{1/2} \to [L^2, H^1]_{1/2} = H^{1/2}, $$

is continuous. In particular, by (4.35), there exists $C > 0$, such that

$$ \forall v \in [W^{-1}, W^0]_{1/2}, \quad \|v\|_{H^{1/2}} \le C \left( \|v\|_{(H_{00}^{1/2}(\text{BL}))'} + \|\nabla v\|_{(H_{00}^{1/2}(\text{BL}))'} \right). \quad (4.39) $$
## 4.6. 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,

$$ \|v\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})} \quad (4.36) $$

satisfied by any distribution $v \in D'(\text{BL})$ such that $v, \nabla v \in H^{-1}(\text{BL})$, the dual of the space $H_0^1(\text{BL})$ as defined by (4.2). Then, we will prove Corollary 4.1, namely

$$ V_1 \subset H^{\frac{1}{2}}(\text{BL}), \quad (4.37) $$

with continuous injection, where $V_1$ is defined by (4.4), equipped with the norm

$$ \|u\|_{V^1} = \|u\|_{L^2(\text{BL})} + \|\nabla_h u\|_{L^2(\text{BL})} + \|z^{1/2} \partial_z u\|_{L^2(\text{BL})}. $$

In particular, being $V_1$ a closed subspace of $H^{1/2}(\text{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 $v \in D'(\text{BL})$ such that $v, \nabla v \in H^{-1}$, and $v_n, \nabla v_n \in V_n^{-1}$ be given by

$$ \forall u \in V_n, U \in V_n^3, \quad \langle v_n, u \rangle = \langle v, u \rangle, \quad \langle \nabla v_n, U \rangle = \langle \nabla v, U \rangle, $$

where $V_n^{-1}$ is defined by (4.11). We have, for all $n \in \mathbb{N}$,

$$ \|\nabla v_n\|_{H^{-1}(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})}, \quad \|v_n\|_{H^{-1}(\text{BL})} \le \|v\|_{H^{-1}(\text{BL})}, $$

and note that (by explicit computation) it easily follows that the constants on the right-hand side are equal to 1, independently of $n \in \mathbb{N}$. Therefore, according to inequality (4.13),

$$ \|v_n\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})}. \quad (4.38) $$

Then, the sequence $(v_n)_{n \in \mathbb{N}}$ is bounded in $L^2(\text{BL})$, and we can extract a subsequence (still denoted by $(v_n)_{n \in \mathbb{N}}$) that weakly converges to some $\tilde{v}$, which is equal to $v$ in $H^{-1}(\text{BL})$, therefore in $L^2(\text{BL})$, which yields $v \in L^2(\text{BL})$ and

$$ \|v\|_{L^2(\text{BL})} \le \liminf_{n \to \infty} \|v_n\|_{L^2(\text{BL})} \le \|\nabla v\|_{H^{-1}(\text{BL})} + \|v\|_{H^{-1}(\text{BL})}, $$

hence (4.36).

*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

$$ \mathrm{Id} : W^{-1} \to L^2 $$

is continuous. Moreover, $\mathrm{Id} : W^0 \to H^1$ is also continuous. Therefore, by the interpolation theorem [12], also the restriction of the identity (denoted still by $\mathrm{Id}$)

$$ \mathrm{Id} : [W^{-1}, W^0]_{1/2} \to [L^2, H^1]_{1/2} = H^{1/2}, $$

is continuous. In particular, by (4.35), there exists $C > 0$, such that

$$ \forall v \in [W^{-1}, W^0]_{1/2}, \quad \|v\|_{H^{1/2}} \le C \left( \|v\|_{(H_{00}^{1/2}(\text{BL}))'} + \|\nabla v\|_{(H_{00}^{1/2}(\text{BL}))'} \right). \quad (4.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

$$V^1 \subset [W^{-1}, W^0]_{1/2}, \qquad (4.40)$$

with continuous injection. In the case relevant for our problem, we have

$$\|z^{1/2}\varphi\|_{L^2} \le C\|\varphi\|_{H_{00}^{1/2}}.$$

Let $u \in V_1$, and $\varphi \in H_{00}^{1/2}(\text{BL})$. Then,

$$\left|\int_{\text{BL}} \partial_z u \varphi \, dx\right| = \left|\int_{\text{BL}} \sqrt{z} \partial_z u \frac{\varphi}{\sqrt{z}} \, dx\right| \le \|\sqrt{z} \partial_z u\|_{L^2} \|z^{-1/2} \varphi\|_{L^2} \le C \|\sqrt{z} \partial_z u\|_{L^2} \|\varphi\|_{H_{00}^{1/2}},$$

hence $\partial_z u \in (H_{00}^{1/2}(\text{BL}))'$, and

$$\|\partial_z u\|_{(H_{00}^{1/2}(\text{BL}))'} \le C \|\sqrt{z} \partial_z u\|_{L^2}.$$

In the same way, we have also

$$\|u\|_{(H_{00}^{1/2}(\text{BL}))'} \le C \sqrt{z_{\text{top}}} \|u\|_{L^2}, \quad \|\nabla_h u\|_{(H_{00}^{1/2}(\text{BL}))'} \le C \sqrt{z_{\text{top}}} \|\nabla_h u\|_{L^2}.$$

Combining all the previous inequalities yields

$$\|u\|_{[W^{-1}, W^0]_{1/2}} \le C\|u\|_{V_1},$$

hence (4.40), which concludes the proof.

<figure><img src="image_2.png" /></figure>

The functional setting for the proof of Theorem 4.1 is more sophisticated than that used in the non-limiting cases $0 \le \alpha < 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 $\alpha = 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 $\alpha = 1$ and a counter-examples is easily built by means of a double logarithmic function. This is why we resort to the space $V_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 $z = 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 $\alpha = 1$. In this case once we have the adapted functional setting we can apply Lax–Milgram in the space $V_1$ with the weak formulation defined by: $\forall v \in V_1$,

$$\lambda \int_{\text{BL}} uv + \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z \partial_z u \partial_z v + C_D \int_{\Gamma_{\text{top}}} uv - \int_{\Gamma_{\text{top}}} v G = \int_{\text{BL}} fv, \quad (4.41)$$

and eventually observe that the proof can be made fully rigorous again by approximation obtained adding $-\epsilon \partial_z^2 u$ to the equations.

<figure><img src="image_3.png" /></figure>

**Remark 4.1.** Note that the trace of $v \in V_1$ is not defined at $z = 0$, where the weight vanishes. On the other hand $v|_{z=z_{\text{top}}}$ is well defined in $H^{1/2}(\Gamma_{\text{top}})$ since the function $v$ belongs to $W^{1,2}$ in a neighborhood of the top part of the boundary. Hence, the integrals $\int_{\Gamma_{\text{top}}} v G$ and $\int_{\Gamma_{\text{top}}} uv$ are properly defined.
*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

$$V^1 \subset [W^{-1}, W^0]_{1/2}, \qquad (4.40)$$

with continuous injection. In the case relevant for our problem, we have

$$\|z^{1/2}\varphi\|_{L^2} \leq C\|\varphi\|_{H_{00}^{1/2}}.$$

Let $u \in V_1$, and $\varphi \in H_{00}^{1/2}(\text{BL})$. Then,

$$\left|\int_{\text{BL}} \partial_z u \varphi \, dx\right| = \left|\int_{\text{BL}} \sqrt{z} \partial_z u \frac{\varphi}{\sqrt{z}} \, dx\right| \leq \|\sqrt{z} \partial_z u\|_{L^2} \|z^{-1/2} \varphi\|_{L^2} \leq C \|\sqrt{z} \partial_z u\|_{L^2} \|\varphi\|_{H_{00}^{1/2}},$$

hence $\partial_z u \in (H_{00}^{1/2}(\text{BL}))'$, and

$$\|\partial_z u\|_{(H_{00}^{1/2}(\text{BL}))'} \leq C \|\sqrt{z} \partial_z u\|_{L^2}.$$

In the same way, we have also

$$\|u\|_{(H_{00}^{1/2}(\text{BL}))'} \leq C \sqrt{z_{\text{top}}} \|u\|_{L^2}, \quad \|\nabla_h u\|_{(H_{00}^{1/2}(\text{BL}))'} \leq C \sqrt{z_{\text{top}}} \|\nabla_h u\|_{L^2}.$$

Combining all the previous inequalities yields

$$\|u\|_{[W^{-1}, W^0]_{1/2}} \leq C \|u\|_{V_1},$$

hence (4.40), which concludes the proof.

<figure><img src="image_3.png" /></figure>

The functional setting for the proof of Theorem 4.1 is more sophisticated than that used in the non-limiting cases $0 \leq \alpha < 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 $\alpha = 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 $\alpha = 1$ and a counter-examples is easily built by means of a double logarithmic function. This is why we resort to the space $V_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 $z = 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 $\alpha = 1$. In this case once we have the adapted functional setting we can apply Lax–Milgram in the space $V_1$ with the weak formulation defined by: $\forall v \in V_1$,

$$\lambda \int_{\text{BL}} uv + \nu_h \int_{\text{BL}} \nabla_h u \cdot \nabla_h v + \mu \int_{\text{BL}} z \partial_z u \partial_z v + C_D \int_{\Gamma_{\text{top}}} uv - \int_{\Gamma_{\text{top}}} v G = \int_{\text{BL}} fv, \quad (4.41)$$

and eventually observe that the proof can be made fully rigorous again by approximation obtained adding $-\epsilon \partial_z^2 u$ to the equations.

<figure><img src="image_4.png" /></figure>

**Remark 4.1.** Note that the trace of $v \in V_1$ is not defined at $z = 0$, where the weight vanishes. On the other hand $v|_{z=z_{\text{top}}}$ is well defined in $H^{1/2}(\Gamma_{\text{top}})$ since the function $v$ belongs to $W^{1,2}$ in a neighborhood of the top part of the boundary. Hence, the integrals $\int_{\Gamma_{\text{top}}} v G$ and $\int_{\Gamma_{\text{top}}} uv$ are properly defined.
Moreover, as is standard for these problems the function $u \in V_1$ can be taken as test function, proving the energy equality

$$ \lambda \|u\|^2 + \nu_h \|\nabla_h u\|^2 + \mu \|\sqrt{z} \partial_z u\|^2 + C_D \int_{\Gamma_{\text{top}}} |u|^2 - \int_{\Gamma_{\text{top}}} G u = \int_{\text{BL}} fu. $$

Moreover, this can be used also to show 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 $\{z = 0\} \cup \{z = z_{\text{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

$$ H_{0,\tau}^1(\text{BL}) = \{u \in H^1(\text{BL}); u = 0 \text{ on } \Gamma_b \text{ and } \partial_z u = 0 \text{ on } \Gamma_{\text{top}}\}, \quad (4.42) $$

instead of $H_0^1(\text{BL})$ as in (4.2), we can easily convert to the case of a Dirichlet problem in a doubled domain $\Pi \times (0, 2z_{\text{top}})$, with a reflection along the line $z = z_{\text{top}}$. Then, the proof remains the same as in the previous case. Note that since we have $V_1$ functions for which the trace is not well-defined, we cannot expect that our weak solution satisfies $u = 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.

## 5. 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 $\alpha$. We solve the problem (1.4) in two dimensions, using the software Freefem++.
We consider a rectangular box $[0, L] \times [0, z_{\text{top}}]$ with $z_{\text{top}} = 100$, $L = 1000$ and we add periodic conditions on the left and right sides $\{0\} \times [0, z_{\text{top}}]$ and $\{L\} \times [0, z_{\text{top}}]$.

The other boundary conditions are Dirichlet at the bottom, and a linear Navier condition at the top involving the roughness coefficient $C_D$:

$$ \begin{cases} \kappa u^* z_{\text{top}} \frac{du}{dz}(x, z_{\text{top}}) = C_D(V(x) - u(x, z_{\text{top}})), & \text{at } z = z_{\text{top}}, \\ u(x, 0) = 0, & \text{at } z = 0, \end{cases} \quad (5.1) $$

where

$$ V(x) = u_{\text{Log}}(z_{\text{top}})(1 + \epsilon(x)), \quad (5.2) $$

and

$$ u_{\text{Log}}(z) = \begin{cases} \frac{C_* u^*}{C_\nu} z, & z \in [0, z_0] \\ \frac{C_* u^*}{C_\nu} \left( \log\left(\frac{z}{z_0}\right) + 1 \right), & z \in [z_0, z_{\text{top}}]. \end{cases} \quad (5.3) $$

The aim of this numerical study is to check if the numerical solution $u$ 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 $u_{\text{Log}}$.

Then, we will discuss the influence of the different parameters ($\alpha$ in Sect. 5.2), $C_D$ and $u^*$ in the Section 5.3 and calculate the difference between the numerical solution and $u_{\text{Log}}$, which allow us in Section 5.4 to derive an analytical formula of the stabilization function $\Psi$ deduced from the numerical results by interpolation.

Finally, we will consider the influence on the $x$-axis of a small perturbations $\epsilon$, as involved in the wind at the top given by (5.2).
Moreover, as is standard for these problems the function $u \in V_1$ can be taken as test function, proving the energy equality

$$ \lambda \|u\|^2 + \nu_h \|\nabla_h u\|^2 + \mu \|\sqrt{z}\partial_z u\|^2 + C_D \int_{\Gamma_{\text{top}}} |u|^2 - \int_{\Gamma_{\text{top}}} G u = \int_{\text{BL}} fu. $$

Moreover, this can be used also to show 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 $\{z = 0\} \cup \{z = z_{\text{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

$$ H_{0,\tau}^1(\text{BL}) = \{u \in H^1(\text{BL}); u = 0 \text{ on } \Gamma_b \text{ and } \partial_z u = 0 \text{ on } \Gamma_{\text{top}}\}, \quad (4.42) $$

instead of $H_0^1(\text{BL})$ as in (4.2), we can easily convert to the case of a Dirichlet problem in a doubled domain $\Pi \times (0, 2z_{\text{top}})$, with a reflection along the line $z = z_{\text{top}}$. Then, the proof remains the same as in the previous case. Note that since we have $V_1$ functions for which the trace is not well-defined, we cannot expect that our weak solution satisfies $u = 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.

## 5. 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 $\alpha$. We solve the problem (1.4) in two dimensions, using the software Freefem++.
We consider a rectangular box $[0, L] \times [0, z_{\text{top}}]$ with $z_{\text{top}} = 100$, $L = 1000$ and we add periodic conditions on the left and right sides $\{0\} \times [0, z_{\text{top}}]$ and $\{L\} \times [0, z_{\text{top}}]$.

The other boundary conditions are Dirichlet at the bottom, and a linear Navier condition at the top involving the roughness coefficient $C_D$:

$$ \begin{cases} \kappa u^* z_{\text{top}} \frac{du}{dz}(x, z_{\text{top}}) = C_D(V(x) - u(x, z_{\text{top}})), & \text{at } z = z_{\text{top}}, \\ u(x, 0) = 0, & \text{at } z = 0, \end{cases} \quad (5.1) $$

where

$$ V(x) = u_{\text{Log}}(z_{\text{top}})(1 + \epsilon(x)), \quad (5.2) $$

and

$$ u_{\text{Log}}(z) = \begin{cases} \frac{C_* u^*}{C_\nu} z, & z \in [0, z_0] \\ \frac{C_* u^*}{C_\nu} \left( \log\left(\frac{z}{z_0}\right) + 1 \right), & z \in [z_0, z_{\text{top}}]. \end{cases} \quad (5.3) $$

The aim of this numerical study is to check if the numerical solution $u$ 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 $u_{\text{Log}}$.

Then, we will discuss the influence of the different parameters ($\alpha$ in Sect. 5.2), $C_D$ and $u^*$ in the Section 5.3 and calculate the difference between the numerical solution and $u_{\text{Log}}$, which allow us in Section 5.4 to derive an analytical formula of the stabilization function $\Psi$ deduced from the numerical results by interpolation.

Finally, we will consider the influence on the $x$-axis of a small perturbations $\epsilon$, as involved in the wind at the top given by (5.2).
<figure><img src="image_4.png" /><figcaption>FIGURE 1. log-law. (a) $z_0 = 0.1$. (b) $z_0 = 1$. (c) $z_0 = 10$.</figcaption></figure>

**Remark 5.1.** We use in the code the command $u = 0$ at $z = 0$, and we take $\lambda = 0$. Contrary to what the analysis predicts, the case $\alpha = 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 $\alpha = 1$. However, we did not have studied yet this numerical aspect of the problem.

## 5.1. Settings of the analysis

### 5.1.1. 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, L] \times [0, z_{\text{top}}]$, where

$$L = 1000 \text{ m and } z_{\text{top}} = 100 \text{ m.}$$

- After several simulations to get velocities that can be measured in situ, it looks like that the best values for the calibration constants $C^*$ and $C_\nu$ were $C^* = 10$ and $C_\nu = 15$, as already mentioned in Section 2.1. Examples of log profiles have been plotted in Figure 1 with these values, $u^* = 10$, $z_{\text{top}} = 100$, and respectively $z_0 = 1$ and $z_0 = 10$.

- The height of the viscous sub-layer $z_0$ in the following simulations will be very small compared to $z_{\text{top}}$, with a ratio smaller than 0.01. Besides, the viscous sub-layer height $z_0$ will taken equal to 0.1, which respects a ratio

$$\frac{z_0}{z_{\text{top}}} = 10^{-3} < 0.01.$$

- The parameter $\lambda$ is chosen equal to 0.

- The source $f$ will be considered as constant: $f = 5$. According to formula (2.7), this means that $\delta T \approx -23$°C for $\beta = 2 \times 10^{-3}$, for instance in the case where the ground is at 0°C and the air is dry and cold at a constant temperature equal to -23°C, a situation that can happen in the mountains.

- The perturbation $\epsilon$ which appears in the expression (5.2) of $V$ will be taken equal to 0 in the next subsections, except in Section 5.5 (see (5.9)).

The velocity constant $u_*$ can be seen as a “wind regime switch” belonging to the speed range $[2 \text{ ms}^{-1}, 10 \text{ 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 $C_D$ will be chosen according to the $u^*$. In the Section 5.3 we will show that $C_D \simeq 10^6$ is giving convincing results.

### 5.1.2. Errors

Let $u_k$ denotes the vertical velocities at $x_k = kL$, given by, for every $z \in [0, z_{\text{top}}]$:

$$u_k(z) = u(x_k, z), \qquad (5.4)$$
<figure><img src="image_5.png" /><figcaption>FIGURE 1. log-law. (a) $z_0 = 0.1$. (b) $z_0 = 1$. (c) $z_0 = 10$.</figcaption></figure>

**Remark 5.1.** We use in the code the command $u = 0$ at $z = 0$, and we take $\lambda = 0$. Contrary to what the analysis predicts, the case $\alpha = 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 $\alpha = 1$. However, we did not have studied yet this numerical aspect of the problem.

## 5.1. Settings of the analysis

### 5.1.1. 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, L] \times [0, z_{\text{top}}]$, where

$$L = 1000 \text{ m and } z_{\text{top}} = 100 \text{ m.}$$

- After several simulations to get velocities that can be measured in situ, it looks like that the best values for the calibration constants $C^*$ and $C_\nu$ were $C^* = 10$ and $C_\nu = 15$, as already mentioned in Section 2.1. Examples of log profiles have been plotted in Figure 1 with these values, $u^* = 10$, $z_{\text{top}} = 100$, and respectively $z_0 = 1$ and $z_0 = 10$.

- The height of the viscous sub-layer $z_0$ in the following simulations will be very small compared to $z_{\text{top}}$, with a ratio smaller than 0.01. Besides, the viscous sub-layer height $z_0$ will taken equal to 0.1, which respects a ratio

$$\frac{z_0}{z_{\text{top}}} = 10^{-3} < 0.01.$$

- The parameter $\lambda$ is chosen equal to 0.

- The source $f$ will be considered as constant: $f = 5$. According to formula (2.7), this means that $\delta T \approx -23^\circ C$ for $\beta = 2 \times 10^{-3}$, for instance in the case where the ground is at $0^\circ C$ and the air is dry and cold at a constant temperature equal to $-23^\circ C$, a situation that can happen in the mountains.

- The perturbation $\epsilon$ which appears in the expression (5.2) of $V$ will be taken equal to 0 in the next subsections, except in Section 5.5 (see (5.9)).

The velocity constant $u_*$ can be seen as a “wind regime switch” belonging to the speed range $[2 \text{ ms}^{-1}, 10 \text{ 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 $C_D$ will be chosen according to the $u^*$. In the Section 5.3 we will show that $C_D \simeq 10^6$ is giving convincing results.

### 5.1.2. Errors

Let $u_k$ denotes the vertical velocities at $x_k = kL$, given by, for every $z \in [0, z_{\text{top}}]$:

$$u_k(z) = u(x_k, z), \qquad (5.4)$$
<figure><img src="image_5.png" /><figcaption>FIGURE 2. Velocity profiles for different horizontal values.</figcaption></figure>

where $k \in \{0, \frac{1}{10}, \frac{2}{10}, \dots, 1\}$. These $u_k$ will be compared with the log profile $u_{\text{Log}}$ defined by (5.3). To achieve these comparisons, we introduce the errors $err$ at a point $x$

$$err(x) = \frac{1}{N} \sum_{i=1}^{N} |u(x, z_i) - u_{\text{Log}}(z_i)|, \qquad (5.5)$$

and the relative error $relerr$ at $x$

$$relerr(x) = \sum_{i=1}^{N} \frac{|u(x, z_i) - u_{\text{Log}}(z_i)|}{|u(x, z_i)|}, \qquad (5.6)$$

which will be more relevant than the error because of the importance of the velocities, when $u^* > 7 \text{ms}^{-1}$ for instance.

Without the perturbation $\epsilon$, the vertical velocities $u_k$ 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

$$\bar{u}(z) := \frac{1}{p} \sum_{k=1}^{p} u(x_k, z), \qquad (5.7)$$

and we will use it instead of the $u_k$. To perform the simulations, $p = 11$ and $x_k$ will belongs to $\{0, 100, 200, \dots, 1000\}$.
<figure><img src="image_6.png" /><figcaption>FIGURE 2. Velocity profiles for different horizontal values.</figcaption></figure>

where $k \in \{0, \frac{1}{10}, \frac{2}{10}, \dots, 1\}$. These $u_k$ will be compared with the log profile $u_{\text{Log}}$ defined by (5.3). To achieve these comparisons, we introduce the errors $err$ at a point $x$

$$err(x) = \frac{1}{N} \sum_{i=1}^{N} |u(x, z_i) - u_{\text{Log}}(z_i)|, \qquad (5.5)$$

and the relative error $relerr$ at $x$

$$relerr(x) = \sum_{i=1}^{N} \frac{|u(x, z_i) - u_{\text{Log}}(z_i)|}{|u(x, z_i)|}, \qquad (5.6)$$

which will be more relevant than the error because of the importance of the velocities, when $u^* > 7 \text{ms}^{-1}$ for instance.

Without the perturbation $\epsilon$, the vertical velocities $u_k$ 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

$$\bar{u}(z) := \frac{1}{p} \sum_{k=1}^{p} u(x_k, z), \qquad (5.7)$$

and we will use it instead of the $u_k$. To perform the simulations, $p = 11$ and $x_k$ will belongs to $\{0, 100, 200, \dots, 1000\}$.
<table><caption>TABLE 1. Errors and relative errors between the vertical velocities $u_k$ and $u_{\text{Log}}$. The parameters taken are $C_D = 10^6$, $\alpha = 1$, $u^* = 4$, $x_1 = 200$, $x_2 = 400$, $x_3 = 600$, $x_4 = 800$.</caption><thead><tr><th colspan="2">u<sub>1</sub></th><th colspan="2">u<sub>2</sub></th><th colspan="2">u<sub>3</sub></th><th colspan="2">u<sub>4</sub></th></tr><tr><th>err</th><th>relerr</th><th>err</th><th>relerr</th><th>err</th><th>relerr</th><th>err</th><th>relerr</th></tr></thead><tbody><tr><td>0.37</td><td>0.042</td><td>0.49</td><td>0.053</td><td>0.53</td><td>0.056</td><td>0.47</td><td>0.047</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>FIGURE 3. Vertical velocities. (a) $\alpha = 0.3$. (b) $\alpha = 0.6$. (c) $\alpha = 0.8$. (d) $\alpha = 1$.</figcaption></figure>

## 5.2. Influence of $\alpha$

We observe that the more $\alpha$ 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 $\bar{u}(z)$, where $z \in \{1, 2, \dots, 100\}$.

As a result, the model is relevant only for $\alpha = 1$. We will take $\alpha = 1$ in the next simulations.

## 5.3. The three different regimes: influence of $C_D$

We have calibrated the constant values to get wind velocities which are physically relevant (in $\text{ms}^{-1}$). We have observed three different regimes for $u_*$ respectively equal to $4 \text{ms}^{-1}$, $7 \text{ms}^{-1}$ and $10 \text{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 $C_D$ coefficients.
<table><caption>TABLE 1. Errors and relative errors between the vertical velocities $u_k$ and $u_{\text{Log}}$. The parameters taken are $C_D = 10^6$, $\alpha = 1$, $u^* = 4$, $x_1 = 200$, $x_2 = 400$, $x_3 = 600$, $x_4 = 800$.</caption><thead><tr><th colspan="2">u<sub>1</sub></th><th colspan="2">u<sub>2</sub></th><th colspan="2">u<sub>3</sub></th><th colspan="2">u<sub>4</sub></th></tr><tr><th>err</th><th>relerr</th><th>err</th><th>relerr</th><th>err</th><th>relerr</th><th>err</th><th>relerr</th></tr></thead><tbody><tr><td>0.37</td><td>0.042</td><td>0.49</td><td>0.053</td><td>0.53</td><td>0.056</td><td>0.47</td><td>0.047</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIGURE 3. Vertical velocities. (a) $\alpha = 0.3$. (b) $\alpha = 0.6$. (c) $\alpha = 0.8$. (d) $\alpha = 1$.</figcaption></figure>

## 5.2. Influence of $\alpha$

We observe that the more $\alpha$ 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 $\bar{u}(z)$, where $z \in \{1, 2, \dots, 100\}$.

As a result, the model is relevant only for $\alpha = 1$. We will take $\alpha = 1$ in the next simulations.

## 5.3. The three different regimes: influence of $C_D$

We have calibrated the constant values to get wind velocities which are physically relevant (in $\text{ms}^{-1}$). We have observed three different regimes for $u_*$ respectively equal to $4 \text{ms}^{-1}$, $7 \text{ms}^{-1}$ and $10 \text{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 $C_D$ coefficients.
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Influence of $C_D$ for $u^* = 10$. (a) $C_D = 100$. (b) $C_D = 1000$. (c) $C_D = 10^6$.</figcaption></figure>

<table><caption>TABLE 2. CD calibration for the different regimes.</caption><thead><tr><th rowspan="2">C<sub>D</sub></th><th colspan="2">u* = 4</th><th colspan="2">u* = 7</th><th colspan="2">u* = 10</th></tr><tr><th>err</th><th>relerr</th><th>err</th><th>relerr</th><th>err</th><th>relerr</th></tr></thead><tbody><tr><td>10<sup>2</sup></td><td>7.58</td><td>0.73</td><td>19.59</td><td>1.59</td><td>32.21</td><td>2.42</td></tr><tr><td>10<sup>3</sup></td><td>1.79</td><td>0.13</td><td>5.50</td><td>0.24</td><td>10.34</td><td>0.32</td></tr><tr><td>10<sup>4</sup></td><td>0.63</td><td>0.061</td><td>1.68</td><td>0.083</td><td>2.85</td><td>0.094</td></tr><tr><td>10<sup>5</sup></td><td>0.50</td><td>0.054</td><td>1.23</td><td>0.067</td><td>1.90</td><td>0.071</td></tr><tr><td>10<sup>6</sup></td><td>0.48</td><td>0.053</td><td>1.18</td><td>0.066</td><td>1.81</td><td>0.069</td></tr><tr><td>10<sup>9</sup></td><td>0.48</td><td>0.053</td><td>1.18</td><td>0.065</td><td>1.79</td><td>0.069</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>FIGURE 5. Errors and relative errors dependence on $C_D$.</figcaption></figure>

We can see in Figure 4 that the vertical velocity $\bar{u}$ in blue is close to the log-law when $C_D$ is big and far when $C_D$ is small. To quantify this, the errors and relative errors corresponding in the Table 2 show that the bigger $C_D$ 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 $u^* = 4$, $u^* = 7$ and $u^* = 10$. The value $10^6$ seems to be the threshold value for the coefficient $C_D$.
<figure><img src="image_8.png" /><figcaption>FIGURE 4. Influence of $C_D$ for $u^* = 10$. (a) $C_D = 100$. (b) $C_D = 1000$. (c) $C_D = 10^6$.</figcaption></figure>

<table><caption>TABLE 2. CD calibration for the different regimes.</caption><thead><tr><th rowspan="2">C<sub>D</sub></th><th colspan="2">u* = 4</th><th colspan="2">u* = 7</th><th colspan="2">u* = 10</th></tr><tr><th>err</th><th>relerr</th><th>err</th><th>relerr</th><th>err</th><th>relerr</th></tr></thead><tbody><tr><td>10<sup>2</sup></td><td>7.58</td><td>0.73</td><td>19.59</td><td>1.59</td><td>32.21</td><td>2.42</td></tr><tr><td>10<sup>3</sup></td><td>1.79</td><td>0.13</td><td>5.50</td><td>0.24</td><td>10.34</td><td>0.32</td></tr><tr><td>10<sup>4</sup></td><td>0.63</td><td>0.061</td><td>1.68</td><td>0.083</td><td>2.85</td><td>0.094</td></tr><tr><td>10<sup>5</sup></td><td>0.50</td><td>0.054</td><td>1.23</td><td>0.067</td><td>1.90</td><td>0.071</td></tr><tr><td>10<sup>6</sup></td><td>0.48</td><td>0.053</td><td>1.18</td><td>0.066</td><td>1.81</td><td>0.069</td></tr><tr><td>10<sup>9</sup></td><td>0.48</td><td>0.053</td><td>1.18</td><td>0.065</td><td>1.79</td><td>0.069</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>FIGURE 5. Errors and relative errors dependence on $C_D$.</figcaption></figure>

We can see in Figure 4 that the vertical velocity $\bar{u}$ in blue is close to the log-law when $C_D$ is big and far when $C_D$ is small. To quantify this, the errors and relative errors corresponding in the Table 2 show that the bigger $C_D$ 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 $u^* = 4$, $u^* = 7$ and $u^* = 10$. The value $10^6$ seems to be the threshold value for the coefficient $C_D$.
<figure><img src="image_9.png" /><figcaption>FIGURE 6. Stabilization functions for different values of $u^*$.</figcaption></figure>

<table><caption>TABLE 3. Errors and relative errors between the stabilized/unstabilized velocities and the log-law.</caption><thead><tr><th>$u^*$</th><th>$C_D$</th><th>err w/o stab</th><th>err with stab</th><th>relerr w/o stab</th><th>relerr with stab</th></tr></thead><tbody><tr><td>4</td><td>$10^6$</td><td>0.49</td><td>0.30</td><td>0.053</td><td>0.026</td></tr><tr><td>7</td><td>$10^6$</td><td>1.18</td><td>0.35</td><td>0.066</td><td>0.021</td></tr><tr><td>10</td><td>$10^6$</td><td>1.81</td><td>0.47</td><td>0.069</td><td>0.020</td></tr></tbody></table>

## 5.4. Stabilization functions

From our numerical results, by an empirical method of extrapolation by successive approximations, we have found the following stabilization function

$$ \Psi(z) = 2u^*(e^{(z_0-z)} - e^{-2z}) - \frac{u^*}{200}z + 0.4u^* \quad (5.8) $$

so that, if $u$ denotes the numerical result,

$$ u + \psi \approx u_{\text{Log}}. $$

It gives for the different values of $u^*$ the curves shown in Figure 6. The peak we can see corresponds to the the height of the viscous sub-layer $z_0$.

We can compare the errors and relative errors between the raw velocity $u$ and the log-law $u_{\text{Log}}$ on the one hand, and between the stabilized velocity $u + \Psi$ and $u_{\text{Log}}$ on the other hand, for the different regimes given by $u^*$. 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 $z_0$.

## 5.5. With horizontal perturbation

In this part we study the effect of a small oscillation in the horizontal direction, given by

$$ \epsilon(x) = 0.01 \sin\left(11 \frac{2\pi x}{L}\right), \quad (5.9) $$

so that $V(x) = u_{\text{Log}}(z_{\text{top}})(1 + \epsilon(x))$.
<figure><img src="image_10.png" /><figcaption>FIGURE 6. Stabilization functions for different values of $u^*$.</figcaption></figure>

<table><caption>TABLE 3. Errors and relative errors between the stabilized/unstabilized velocities and the log-law.</caption><thead><tr><th>u<sup>*</sup></th><th>C<sub>D</sub></th><th>err w/o stab</th><th>err with stab</th><th>relerr w/o stab</th><th>relerr with stab</th></tr></thead><tbody><tr><td>4</td><td>10<sup>6</sup></td><td>0.49</td><td>0.30</td><td>0.053</td><td>0.026</td></tr><tr><td>7</td><td>10<sup>6</sup></td><td>1.18</td><td>0.35</td><td>0.066</td><td>0.021</td></tr><tr><td>10</td><td>10<sup>6</sup></td><td>1.81</td><td>0.47</td><td>0.069</td><td>0.020</td></tr></tbody></table>

## 5.4. Stabilization functions

From our numerical results, by an empirical method of extrapolation by successive approximations, we have found the following stabilization function

$$ \Psi(z) = 2u^*(e^{(z_0-z)} - e^{-2z}) - \frac{u^*}{200}z + 0.4u^* \quad (5.8) $$

so that, if $u$ denotes the numerical result,

$$ u + \psi \approx u_{\text{Log}}. $$

It gives for the different values of $u^*$ the curves shown in Figure 6. The peak we can see corresponds to the height of the viscous sub-layer $z_0$.

We can compare the errors and relative errors between the raw velocity $u$ and the log-law $u_{\text{Log}}$ on the one hand, and between the stabilized velocity $u + \Psi$ and $u_{\text{Log}}$ on the other hand, for the different regimes given by $u^*$. 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 $z_0$.

## 5.5. With horizontal perturbation

In this part we study the effect of a small oscillation in the horizontal direction, given by

$$ \epsilon(x) = 0.01 \sin\left(11 \frac{2\pi x}{L}\right), \quad (5.9) $$

so that $V(x) = u_{\text{Log}}(z_{\text{top}})(1 + \epsilon(x))$.
<figure><img src="image_10.png" /><figcaption>FIGURE 7. Velocity profile with and without stabilization compared with the log-law. (a) $u^* = 4$. (b) $u^* = 7$. (c) $u^* = 10$.</figcaption></figure>

<table><caption>TABLE 4. Errors and relative errors between the velocities and the log-law.</caption><thead><tr><th rowspan="2">u<sup>*</sup></th><th colspan="3">err</th><th colspan="3">relerr</th></tr><tr><th>w/o &epsilon;</th><th>with &epsilon;</th><th>with &epsilon; and &Psi;</th><th>w/o &epsilon;</th><th>with &epsilon;</th><th>with &epsilon; and &Psi;</th></tr></thead><tbody><tr><td>4</td><td>0.49</td><td>0.62</td><td>0.27</td><td>0.053</td><td>0.061</td><td>0.024</td></tr><tr><td>7</td><td>1.18</td><td>1.42</td><td>0.41</td><td>0.066</td><td>0.073</td><td>0.023</td></tr><tr><td>10</td><td>1.81</td><td>2.15</td><td>0.65</td><td>0.069</td><td>0.078</td><td>0.024</td></tr></tbody></table>

We keep the value $C_D = 10^6$, $\alpha = 1$, and we will compare the vertical velocities $\bar{u}$ we get with the ones we had without this perturbation for the three different $u^*$. 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 $\Psi$, we can manage to keep the errors still small, even with the perturbation $\epsilon$. This opens an interesting stability problem à la “Lyapunov”.

## 5.6. Conclusion and perspectives

We have shown numerically that the solution $u$ of the problem (1.4) is a very good approximation for the known Monin-Obukhov log-law when $\alpha = 1$ for adiabatic flows, with a large calibration constant $C_D$, imposing a Dirichlet condition at $\Gamma_{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 $T$, where the source term $f$ is given by (2.7), which is a work under progress.

### Acknowledgements

L.C Berselli acknowledges support by INdAM-GNAMPA and by MIUR, within the project PRIN20204NT8W4.004 Nonlinear evolution PDEs, fluid dynamics and transport equations: theoretical foundations and applications.

## REFERENCES

[1] C. Amrouche, P.G. Ciarlet and C. Mardare, On a lemma of Jacques-Louis Lions and its relation to other fundamental results. *J. Math. Pures Appl.* **104** (2015) 207–226.

[2] C. Amrouche, L.C. Berselli, R. Lewandowski and D.D. Nguyen, Turbulent flows as generalized Kelvin-Voigt materials: modeling and analysis. *Nonlinear Anal.* **196** (2020) 111790.

[3] L.C. Berselli and D. Breit, On the existence of weak solutions for the steady Baldwin-Lomax model and generalizations. *J. Math. Anal. App.* **501** (2021) 124633.
<figure><img src="image_11.png" /><figcaption>FIGURE 7. Velocity profile with and without stabilization compared with the log-law. (a) $u^* = 4$. (b) $u^* = 7$. (c) $u^* = 10$.</figcaption></figure>

<table><caption>TABLE 4. Errors and relative errors between the velocities and the log-law.</caption><thead><tr><th rowspan="2">u<sup>*</sup></th><th colspan="3">err</th><th colspan="3">relerr</th></tr><tr><th>w/o ε</th><th>with ε</th><th>with ε and Ψ</th><th>w/o ε</th><th>with ε</th><th>with ε and Ψ</th></tr></thead><tbody><tr><td>4</td><td>0.49</td><td>0.62</td><td>0.27</td><td>0.053</td><td>0.061</td><td>0.024</td></tr><tr><td>7</td><td>1.18</td><td>1.42</td><td>0.41</td><td>0.066</td><td>0.073</td><td>0.023</td></tr><tr><td>10</td><td>1.81</td><td>2.15</td><td>0.65</td><td>0.069</td><td>0.078</td><td>0.024</td></tr></tbody></table>

We keep the value $C_D = 10^6$, $\alpha = 1$, and we will compare the vertical velocities $\bar{u}$ we get with the ones we had without this perturbation for the three different $u^*$. 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 $\Psi$, we can manage to keep the errors still small, even with the perturbation $\epsilon$. This opens an interesting stability problem à la “Lyapunov”.

## 5.6. Conclusion and perspectives

We have shown numerically that the solution $u$ of the problem (1.4) is a very good approximation for the known Monin-Obukhov log-law when $\alpha = 1$ for adiabatic flows, with a large calibration constant $C_D$, imposing a Dirichlet condition at $\Gamma_{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 $T$, where the source term $f$ is given by (2.7), which is a work under progress.

### Acknowledgements

L.C Berselli acknowledges support by INdAM-GNAMPA and by MIUR, within the project PRIN20204NT8W4.004 Nonlinear evolution PDEs, fluid dynamics and transport equations: theoretical foundations and applications.

## REFERENCES

[1] C. Amrouche, P.G. Ciarlet and C. Mardare, On a lemma of Jacques-Louis Lions and its relation to other fundamental results. *J. Math. Pures Appl.* **104** (2015) 207–226.

[2] C. Amrouche, L.C. Berselli, R. Lewandowski and D.D. Nguyen, Turbulent flows as generalized Kelvin-Voigt materials: modeling and analysis. *Nonlinear Anal.* **196** (2020) 111790.

[3] L.C. Berselli and D. Breit, On the existence of weak solutions for the steady Baldwin-Lomax model and generalizations. *J. Math. Anal. App.* **501** (2021) 124633.
[4] L.C. Berselli, A. Kaltenbach, R. Lewandowski and M. Růžička, On the existence of weak solutions for a family of unsteady rotational smagorinsky models. To appear in *Pure Appl. Funct. Anal.* arXiv:2107.00236 (2023).

[5] J.W. Deardorff, Numerical investigation of neutral and unstable planetary boundary layers. *J. Atmos. Sci.* **29** (1972) 91–115.

[6] G. Duvaut and J.-L. Lions, Inequalities in Mechanics and Physics. Vol. 219 of *Grundlehren der Mathematischen Wissenschaften*. Springer-Verlag, Berlin-New York (1976). Translated from the French by C. W. John.

[7] F. Hecht, New development in FreeFem++. *J. Numer. Math.* **20** (2012) 251–265.

[8] A. Kufner, Weighted Sobolev Spaces. Vol. 31 of *Teubner-Texte zur Mathematik [Teubner Texts in Mathematics]*. BSB B. G. Teubner Verlagsgesellschaft, Leipzig (1980). With German, French and Russian summaries.

[9] F. Legeais and R. Lewandowski, Continuous boundary condition at the interface for two coupled fluids. *Appl. Math. Lett.* **135** (2023) 108393.

[10] R. Lewandowski, Analyse mathématique et océanographie. Vol. 39 of *Recherches en Mathématiques Appliquées [Research in Applied Mathematics]*. Masson, Paris (1997).

[11] R. Lewandowski and G. Pichot, Numerical simulation of water flow around a rigid fishing net. *Comput. Methods Appl. Mech. Eng.* **196** (2007) 4737–4754.

[12] J.-L. Lions and E. Magenes, Non-Homogeneous Boundary Value Problems and Applications. Vol. I. *Die Grundlehren der mathematischen Wissenschaften, Band 181*. Springer-Verlag, New York-Heidelberg (1972). Translated from the French by P. Kenneth.

[13] B. Mohammadi and O. Pironneau, Analysis of the *k*-epsilon Turbulence Model. RAM: Research in Applied Mathematics. Masson, Paris; John Wiley & Sons, Ltd., Chichester (1994).

[14] A.S. Monin and A.M. Obukhov, Basic laws of turbulent mixing in the surface layer of the atmosphere. *Tr. Akad. Nauk. SSSR Geophiz. Inst.* **24** (1954) 163–187.

[15] A.M. Obukhov, Turbulence in an atmosphere with a non-uniform temperature. *Boundary-Layer Meteorol.* **2** (1971) 7–29.

[16] C. Pelletier, *Mathematical study of the air-sea coupling problem including turbulent scale effects*. Thesis, University of Grenoble, France, UCLouvain, Louvain-la-Neuve, Belgium (2018).

[17] B. Pinier, E. Mémin, S. Laizet and R. Lewandowski, Stochastic flow approach to model the mean velocity profile of wall-bounded flows. *Phys. Rev. E* **99** (2019) 063101.

[18] B. Pinier, R. Lewandowski, E. Mémin and P. Chandramouli, Testing a one-closure equation turbulence model in neutral boundary layers. *Comput. Methods Appl. Mech. Eng.* **376** (2021) 113662.

[19] L. Prandtl, Motion of fluids with very little viscosity, in *Verhandlungen des dritten internationalen Mathematiker-Kongresses in Heidelberg* 1904, p.484, edited by A. Krazer. Teubner, Leipzig, Germany (1905).

[20] J. Rappaz and J. Rochat, On non-linear Stokes problems with viscosity depending on the distance to the wall. *C. R. Math. Acad. Sci. Paris* **354** (2016) 499–502.

[21] J. Rappaz and J. Rochat, On some weighted Stokes problems: applications on Smagorinsky models, in *Contributions to Partial Differential Equations and Applications*. Vol. 47 of *Comput. Methods Appl. Sci.* Springer, Cham (2019) 395–410.

[22] J. Rappaz and J. Rochat, On von Karman modeling for turbulent flow near a wall. *Methods Appl. Anal.* **26** (2019) 291–295.

[23] T.C. Rebollo and R. Lewandowski, A variational finite element model for large-eddy simulations of turbulent flows. *Chinese Ann. Math. Ser. B* **34** (2013) 667–682.

[24] T.C. Rebollo and R. Lewandowski, Mathematical and Numerical Foundations of Turbulence Models and Applications. *Modeling and Simulation in Science, Engineering and Technology*. Birkhäuser/Springer, New York (2014).

[25] H. Schlichting, Boundary Layer Theory. *McGraw-Hill Series in Mechanical Engineering*. McGraw-Hill Book Co., Inc., New York (1979).

[26] T. von Kármán, Mechanische Ähnlichkeit und turbulenz. *Nachr. Ges. Wiss. Göttingen Math. Phys. Klasse.* **58** (1930) 271–286.

Please help to maintain this journal in open access!

This journal is currently published in open access under the *Subscribe to Open model (S2O)*. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.

<figure><img src="image_11.png" /></figure>
[4] L.C. Berselli, A. Kaltenbach, R. Lewandowski and M. Růžička, On the existence of weak solutions for a family of unsteady rotational smagorinsky models. To appear in *Pure Appl. Funct. Anal.* arXiv:2107.00236 (2023).

[5] J.W. Deardorff, Numerical investigation of neutral and unstable planetary boundary layers. *J. Atmos. Sci.* **29** (1972) 91–115.

[6] G. Duvaut and J.-L. Lions, Inequalities in Mechanics and Physics. Vol. 219 of *Grundlehren der Mathematischen Wissenschaften*. Springer-Verlag, Berlin-New York (1976). Translated from the French by C. W. John.

[7] F. Hecht, New development in FreeFem++. *J. Numer. Math.* **20** (2012) 251–265.

[8] A. Kufner, Weighted Sobolev Spaces. Vol. 31 of *Teubner-Texte zur Mathematik [Teubner Texts in Mathematics]*. BSB B. G. Teubner Verlagsgesellschaft, Leipzig (1980). With German, French and Russian summaries.

[9] F. Legeais and R. Lewandowski, Continuous boundary condition at the interface for two coupled fluids. *Appl. Math. Lett.* **135** (2023) 108393.

[10] R. Lewandowski, Analyse mathématique et océanographie. Vol. 39 of *Recherches en Mathématiques Appliquées [Research in Applied Mathematics]*. Masson, Paris (1997).

[11] R. Lewandowski and G. Pichot, Numerical simulation of water flow around a rigid fishing net. *Comput. Methods Appl. Mech. Eng.* **196** (2007) 4737–4754.

[12] J.-L. Lions and E. Magenes, Non-Homogeneous Boundary Value Problems and Applications. Vol. I. *Die Grundlehren der mathematischen Wissenschaften, Band 181*. Springer-Verlag, New York-Heidelberg (1972). Translated from the French by P. Kenneth.

[13] B. Mohammadi and O. Pironneau, Analysis of the *k*-epsilon Turbulence Model. RAM: Research in Applied Mathematics. Masson, Paris; John Wiley & Sons, Ltd., Chichester (1994).

[14] A.S. Monin and A.M. Obukhov, Basic laws of turbulent mixing in the surface layer of the atmosphere. *Tr. Akad. Nauk. SSSR Geophiz. Inst.* **24** (1954) 163–187.

[15] A.M. Obukhov, Turbulence in an atmosphere with a non-uniform temperature. *Boundary-Layer Meteorol.* **2** (1971) 7–29.

[16] C. Pelletier, *Mathematical study of the air-sea coupling problem including turbulent scale effects*. Thesis, University of Grenoble, France, UCLouvain, Louvain-la-Neuve, Belgium (2018).

[17] B. Pinier, E. Mémin, S. Laizet and R. Lewandowski, Stochastic flow approach to model the mean velocity profile of wall-bounded flows. *Phys. Rev. E* **99** (2019) 063101.

[18] B. Pinier, R. Lewandowski, E. Mémin and P. Chandramouli, Testing a one-closure equation turbulence model in neutral boundary layers. *Comput. Methods Appl. Mech. Eng.* **376** (2021) 113662.

[19] L. Prandtl, Motion of fluids with very little viscosity, in *Verhandlungen des dritten internationalen Mathematiker-Kongresses in Heidelberg* 1904, p.484, edited by A. Krazer. Teubner, Leipzig, Germany (1905).

[20] J. Rappaz and J. Rochat, On non-linear Stokes problems with viscosity depending on the distance to the wall. *C. R. Math. Acad. Sci. Paris* **354** (2016) 499–502.

[21] J. Rappaz and J. Rochat, On some weighted Stokes problems: applications on Smagorinsky models, in *Contributions to Partial Differential Equations and Applications*. Vol. 47 of *Comput. Methods Appl. Sci.* Springer, Cham (2019) 395–410.

[22] J. Rappaz and J. Rochat, On von Karman modeling for turbulent flow near a wall. *Methods Appl. Anal.* **26** (2019) 291–295.

[23] T.C. Rebollo and R. Lewandowski, A variational finite element model for large-eddy simulations of turbulent flows. *Chinese Ann. Math. Ser. B* **34** (2013) 667–682.

[24] T.C. Rebollo and R. Lewandowski, Mathematical and Numerical Foundations of Turbulence Models and Applications. *Modeling and Simulation in Science, Engineering and Technology*. Birkhäuser/Springer, New York (2014).

[25] H. Schlichting, Boundary Layer Theory. *McGraw-Hill Series in Mechanical Engineering*. McGraw-Hill Book Co., Inc., New York (1979).

[26] T. von Kármán, Mechanische Ähnlichkeit und turbulenz. *Nachr. Ges. Wiss. Göttingen Math. Phys. Klasse.* **58** (1930) 271–286.

Please help to maintain this journal in open access!

This journal is currently published in open access under the *Subscribe to Open model (S2O)*. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.

<figure><img src="image_12.png" /></figure>