# Towards a finite volume discretization of the atmospheric surface layer consistent with physical theory

Simon Clément, Florian Lemarié, Eric Blayo

► To cite this version:

Simon Clément, Florian Lemarié, Eric Blayo. Towards a finite volume discretization of the atmospheric surface layer consistent with physical theory. FVCA 2023 - Finite Volumes for Complex Applications 10, Oct 2023, Strasbourg, France. pp.1-8, 10.1007/978-3-031-40864-9. hal-04097172

HAL Id: hal-04097172

https://inria.hal.science/hal-04097172

Submitted on 15 May 2023

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
# Towards a finite volume discretization of the atmospheric surface layer consistent with physical theory

Simon Clément¹, Florian Lemarié¹, and Eric Blayo²

¹ Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

² Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000 Grenoble, France

**Abstract.** We study an atmospheric column and its discretization. Because of numerical considerations, the column must be divided into two parts: (1) a surface layer, excluded from the computational domain and parameterized, and (2) the rest of the column, which reacts more slowly to variations in surface conditions. A usual practice in atmospheric models is to parameterize the surface layer without excluding it from the computational domain, leading to possible consistency issues. We propose here to unify the two representations in a Finite Volume discretization. In order to do so, the reconstruction inside the first grid cell is performed using the particular functions involved in the parameterizations and not only with polynomials. Using a consistency criterion, surface layer management strategies are compared in different physical situations.

**Keywords:** Finite Volume, Monin-Obukhov theory, Surface flux scheme

## 1 Introduction

A common difficulty for atmospheric models is to represent the surface layer (SL), i.e. the area directly and almost instantaneously influenced by the presence of the ground or the ocean. The scales in the SL (approximately the first 10 meters of the air column) are so small that the resolution needed for a numerical model to represent the phenomena correctly is out of reach. However, the Monin-Obukhov (MO) theory, which generalises the wall law to density-stratified fluids, provides under certain simple hypotheses (quasi-stationarity, horizontal homogeneity, etc.) an analytical formulation of the solution in the SL and the expression of the fluxes (heat, momentum) exchanged with the atmosphere above it. At the discrete level, the present treatment of this SL in numerical models is inconsistent: it is both treated like the rest of the atmosphere column by a standard numerical scheme (polynomial profile) when discretising the equations, and in a parameterised form (MO profile, which is a perturbation of a logarithmic profile, see e.g. [1]) for the calculation of fluxes. The consequences of this inconsistency are still poorly assessed in the context of atmospheric modeling, but we can mention for instance that in the context of combustion, it is mentioned in [2] that the way the wall law is implemented in a given code and the
way it interacts with the numerical methods used (in particular the turbulence scheme) can influence the numerical results as much as a particular choice of wall law. In this paper, we will address this inconsistency and propose a new finite volume formulation to remedy it.

**The turbulent Ekman layer model** Our approach is derived hereafter in the case of the 1D vertical Ekman layer model [3] in the neutral case. It includes the Coriolis effect with a constant parameter $f$ and a vertical turbulent flux term $\langle w'u' \rangle$:

$$ \partial_t u + ifu + \partial_z \langle w'u' \rangle = ifu_G \quad (1) $$

where $i$ is the imaginary unit. The constant nudging term $u_G$ pulls the solution towards the geostrophic equilibrium (a large-scale solution where the pressure gradient balances the Coriolis force). The horizontal wind $u$ (in m.s$^{-1}$) is a complex variable accounting for both orientation and speed of the wind. The so-called Boussinesq hypothesis states that the turbulent flux is proportional to the gradient of $u$: $\langle w'u' \rangle = -K_u \partial_z u$ where $K_u \ge 0$ is the turbulent viscosity. In the SL, the MO theory states that this turbulent flux is constant along the vertical axis and provides its analytical expression. We thus obtain the system of equations given in Fig. 1.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Continuous equations in the computational domain ($\delta_a, +\infty$) and constant flux in the SL (0, $\delta_a$). The MO theory specifies the complex constant $u_*e_\tau$ and a range of possible values for $\delta_a$.</figcaption></figure>

**Usual approach in atmospheric models.** The usual choice made in current atmospheric models is to consider that the SL extends from the wall to the center of the first cell (let us note $z_{1/2}$ this altitude). In practice this means that compatibility constraints should apply at $z = z_{1/2}$ to correctly connect the profile as parameterized in the SL with the upper profile obtained using the numerical model. However the usual practice is to integrate the region $z \le z_{1/2}$ corresponding to the SL into the computational domain and then use the MO parameterization only to predict a surface flux at $z = 0$. The impact of this approximation is poorly documented to date: (i) the extent of the SL is fixed for purely numerical reasons and not for physical reasons, which does not guarantee that the solution will converge with the resolution; (ii) the solution in the area $z \le z_{1/2}$ is both parameterized and computed by the model, without ensuring the consistency of these two profiles. The coupling between the SL and the rest of the model is thus weak: the model provides the flow information at $z = z_{1/2}$ to the SL scheme and the latter provides in exchange a surface flux to the model at $z = 0$. In general no other interaction exists. For example, with this kind of coupling, the SL
structures cannot really interact with the rest of the flow. For more details and numerous references, see [4].

A few studies address this issue: several alternatives for implementing a wall law in a Large Eddy Simulation solver are proposed in [2]; a first step toward a proper Finite Volume (FV) approach is proposed by [5], where the authors extend the SL to the entire first cell (0, z₁) and design a scheme adapted to FV.

In this paper, we propose to implement directly in the FV discretization the existing assumptions underlying the MO theory. In order to do so, the reconstruction inside the first grid cell is performed using the analytical functions involved in the wall laws and not only with polynomials. This approach also allows to relax the artificial assumption $\delta_a = z_{1/2}$ and to extend the height of the SL beyond the first grid point if necessary (note that the log-layer mismatch, a well known numerical problem in Large Eddy Simulations, comes from a too thin SL). By being able to choose the thickness of the SL based on physical — and not only numerical — criteria, the consistency of the schemes is improved. The vertical resolution can thus be refined without changing the continuous equations solved by the discretization, thus answering the issues raised in [6]. Numerical experiments with a 1D Ekman layer model are performed, and SL management strategies are compared for different types of stratification.

## 2 The Finite Volume scheme

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Summary of the notations related to the discretisation.</figcaption></figure>

*Spline reconstruction of solutions.* The space domain is divided into M cells delimited by heights ($z_0 = 0, \dots, z_m, \dots, z_M$). The size of the m-th cell is $h_{m+\frac{1}{2}} = z_{m+1} - z_m$ and the average of $u(z)$ over this cell is noted $\bar{u}_{m+\frac{1}{2}} = \frac{1}{h_{m+\frac{1}{2}}} \int_{z_m}^{z_{m+1}} u(z) dz$. The space derivative of $u$ at $z_m$ is noted $\phi_m$. Fig. 2 summarizes these notations. Averaging the evolution equation over a cell gives the semi-discrete equation

$$ (\partial_t + i\bar{f})\bar{u}_{m+\frac{1}{2}} - \frac{K_{u,m+1}\phi_{m+1} - K_{u,m}\phi_m}{h_{m+\frac{1}{2}}} = i\bar{f}u_G \quad (2) $$
The reconstruction of $u(z) = S_{m+\frac{1}{2}}(z-z_{m+\frac{1}{2}})$ is chosen to be a quadratic polynomial (higher order schemes can be similarly derived, see [1]). The continuity of $u(z)$ and its space derivative $\phi$ between cells yields the relation:

$$ \frac{h_{m-1/2}}{6} \phi_{m-1} + \frac{h_{m+1/2} + h_{m-1/2}}{3} \phi_m + \frac{h_{m+1/2}}{6} \phi_{m+1} = \bar{u}_{m+\frac{1}{2}} - \bar{u}_{m-\frac{1}{2}} \quad (3) $$

which is a FV approximation used in fourth-order compact schemes for the first derivative $\partial_z u$ and second-order for $\partial_z^2 u$ (e.g. [7]).

*Usual treatment of the SL with Finite Volume methods.* The typical treatment of the SL in atmospheric models is to use the evolution equation in the first cell ($z_0, z_1$) to compute $\bar{u}_{\frac{1}{2}}$ and then assume that this averaged value is the wind speed at the center of the cell in the Monin-Obukhov theory applied with $\delta_a = z_1$. The corresponding bottom boundary condition is then

$$ \underbrace{K_{u,0} \phi_0^{n+1}}_{\text{Surface flux}} = u_*^2 e_{\tau} \quad \text{with } u_* = \text{BULK}(\underbrace{\bar{u}_{\frac{1}{2}}^n}_{\text{Average around } z_{\frac{1}{2}}}) \quad (4) $$

where BULK is a routine based on the Monin-Obukhov theory, $e_{\tau} = \frac{\bar{u}_{\frac{1}{2}}^{n+1}}{\|\bar{u}_{\frac{1}{2}}^n\|}$, and $n$ denotes the time step. This method has several drawbacks:

- The value at the center of the cell is systematically larger than the average value because of the concavity of Monin-Obukhov profiles. This leads to a systematic underestimation of the surface flux by the SL scheme. A specific SL scheme was designed in [5] to prevent this bias.

- The evolution equation is not compatible with the constant flux hypothesis that defines the SL, as introduced in Section 1.

- $\delta_a$, the height of the SL, is driven only by the space step and does not take into account any physical consideration.

*On the incompatibility.* According to the wall law, $K_{u,0}$ should be equal to the (very small) molecular viscosity $K_{mol} \approx 10^{-5} \text{ m}^2 \cdot \text{s}^{-1}$. However, the boundary condition $K_{u,0}\phi_0 = u_*^2 e_{\tau}$ does not have the same influence depending on the numerical scheme used to discretize (2):

- **Finite Differences:** Injecting the boundary condition in the evolution equation at the first grid level gives

$$ (\partial_t + i\bar{f})u_{1/2} = \frac{1}{h_{1/2}} \left( K_{u,1} \frac{u_{3/2} - u_{1/2}}{h_1} - u_*^2 e_{\tau} \right) \quad (5) $$

where one can see that the value $K_{u,0}$ does not intervene in the equation.

- **Finite Volumes:** applying $(\partial_t + i\bar{f})$ to (3) and using the polynomial reconstruction and the equations of Fig. 1, one can see that the FV scheme implicitly uses

$$ (\partial_t + i\bar{f})u(z_1) = \frac{K_{u,1}\phi_1 - u_*^2 e_{\tau}}{h_{1/2}} + (\partial_t + i\bar{f}) \left( \frac{\phi_1}{3} + \frac{u_*^2 e_{\tau}}{6K_{u,0}} \right) h_{1/2} \quad (6) $$
The (small) value of $K_{u,0}$ directly appears when we assume the parabolic profile inside the first grid cell. As a result, $u(z_1)$ scales with $\frac{1}{K_0}$ and exhibits unreasonable values. To obtain physically plausible profiles, one can replace $K_{u,0}$ by $K_{u,\delta}$: the wall law is then denied and $(\partial_z u)(z_0)$ is multiplied by $\frac{K_{mol}}{K_{u,\delta}}$. Note that this problem would not occur if the simple FV approximation $h_m \phi_m \approx \bar{u}_{m+\frac{1}{2}} - \bar{u}_{m-\frac{1}{2}}$ was used instead of (3).

*Toward a Finite Volume scheme coherent with the physical theory.* To address the drawbacks of the usual method presented above, we now construct a numerical boundary condition that is coherent with the continuous model with a free value of $\delta_a$, named “FV free”:

$$ \underbrace{K_{u,\delta} \phi_{\delta}^{n+1}}_{\text{Flux at } \delta_a} = u_{\star}^2 e_{\tau}^{\text{free}} \quad \text{with } u_{\star} = \text{BULK}(\underbrace{u^n(\delta_a)}_{\text{Reconstruction at } \delta_a}) \quad (7) $$

where $e_{\tau}^{\text{free}} = \frac{u^{n+1}(\delta_a)}{\|u^n(\delta_a)\|}$ is the orientation of $u(\delta_a)$ obtained with the spline reconstruction. For the sake of simplicity, we assume in the following that $\delta_a < z_1$ (this hypothesis being easily relaxed by using the Monin-Obukhov profiles as the reconstruction in cells entirely contained in the SL). In the first grid cell, we assume that the constant flux hypothesis applies for $z < \delta_a$ and we separate this cell into two parts: the surface layer $(0, \delta_a)$ and the “sub-cell” $(\delta_a, z_1)$. This split corresponds to the change of governing equations in Fig. 1. Let $\tilde{h} = z_1 - \delta_a$ be the size of the upper sub-cell $(\delta_a, z_1)$ and $\tilde{u} = \frac{1}{\tilde{h}} \int_{\delta_a}^{z_1} u(z) dz$ be the corresponding averaged value of $u$. The following subgrid reconstruction is used:

$$ u(z) = \begin{cases} S_{1/2} \left( z - \frac{z_1 + \delta_a}{2} \right), & z \ge \delta_a \\ \int_0^z \frac{u_{\star}^2 e_{\tau}^{\text{free}}}{K_{u,z'}} dz', & z < \delta_a \end{cases} \quad (8) $$

where a closed-form of the integral for $z < \delta_a$ is given by MO theory. The quadratic spline $S_{1/2}$ used for reconstruction is computed with the averaged value $\tilde{u}$, the size of the sub-cell $\tilde{h}$ and the fluxes at the extremities $\phi_\delta$ and $\phi_1$: its definition $S_{1/2}(\xi) = \tilde{u} + \frac{\phi_1+\phi_\delta}{2}\xi + \frac{\phi_1-\phi_\delta}{2\tilde{h}}\left(\xi^2 - \frac{\tilde{h}^2}{12}\right)$ is thus similar to the one in the other cells.

## 3 Numerical Experiments

The strategies to handle the SL are now compared through a test of consistency: for several strategies, the differences between a low-resolution and a high-resolution simulations are compared. The smaller the difference between the low-resolution and the high-resolution simulations, the better the consistency of the scheme. The proposed strategy “FV free” is compared with “FV1” (the typical current practice with Finite Volumes), “FV2” (an intermediate between “FV free” and “FV1”: similar to “FV free” but where the height of the surface layer $\delta_a$ is set to $z_1$) and “FD” (a Finite Difference reference).
- The turbulent viscosity is parameterized with a one-equation turbulence closure based on turbulent kinetic energy. The code is available at [8] and an in-depth description in [1]. An Euler implicit time scheme integrates the model over a full day of simulation.

- Parameters are $\Delta t = 30$ s, $u_G = 8$ m.s$^{-1}$, $f = 10^{-4}$ s$^{-1}$

- For “FV free” the same $\delta_a$ is used in both low- and high-resolution simulations, whereas the resolution imposes $\delta_a$ in the other configurations.

- The vertical levels of the low resolution simulation are taken as the 25 first of the 137-level configuration of the atmospheric model *Integrated Forecasting System* at ECMWF (European Centre for Medium-Range Weather Forecasts). The high-resolution simulation has 3 times more cells: two grid levels are added in each of the low-resolution grid cells. The usual SL strategies are designed for low-resolution configurations: the latter can hence be considered as reference solutions, compared through the sensitivity to the resolution.

**Neutral case:** In the neutral case (constant density), the difference between

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Relative difference between low-resolution and high-resolution simulations with several strategies for the SL handling. Left: Relative difference in $u_\star$ as a function of time. Center: vertical profiles of the wind speed at the end of the simulation. Right: Relative difference of the wind speed between low- and high-resolution along the vertical (note the log scale).</figcaption></figure>

low and high resolution of the “FV free” scheme is small at low altitude (see Fig. 3). This is mainly due to two factors:

- $\delta_a = z_{\frac{1}{2} \text{low-res}}$ is the same for both low and high resolutions, whereas for the other surface flux schemes the continuous equations change with $\delta_a$.
- The initial relative difference for $u_*$ (Fig. 3, left panel) is already much smaller than with the other schemes. This is a consequence of the imposed wall law: at initialization, there is already a logarithmic profile in the surface layer, instead of evolving toward a kind of compromise between the parameterized and the modeled values.

**Stratified case:** We now focus on a stratified model [3] that includes more of the physical behavior of atmospheric models: the turbulence closure depends on the density $\rho$ such that $\partial_z\rho \propto -\partial_z\theta$ where $\theta$ is the potential temperature. We designed two cases:

1. A stable stratification, obtained with an initial temperature increasing with the altitude, and a surface temperature decreasing with time. The initial potential temperature is 265 K in the first 100 meters of the atmosphere and then gains 1 degree every 100 meters; the surface temperature starts at 265 K and loses 1 degree every ten hours. The “low resolution” uses 15 grid points in the 400m column and the “high resolution” uses 45 grid points.

2. An unstable stratification, obtained with a surface temperature following a daily oscillation between 279 K and 281 K, and initial profiles of temperature and wind set to constant values of 280 K and $8 \text{ m.s}^{-1}$ respectively. The “low resolution” is composed of 50 grid levels of 10m each; 15 additional stretched levels between 500m and 1080m make sure that the upper boundary condition is not involved in the results. The “high resolution” divides every space cell into 3 new space cells of equal sizes.

The differences between the two simulations are displayed in Fig. 4. In the stable case, the difference between the high resolution and the low resolution results does not significantly change with the surface flux schemes. The “FV2” scheme is not very consistent because it tries to follow the continuous model but with $\delta_a$ changing with the resolution. The Finite Difference or the “FV1” methods suffer less from this problem because, even if $\delta_a$ changes, it is assumed that the evolution equation is integrated inside the surface layer. In [9], authors also find that the sensitivity of their Large Eddy Simulation model to the grid spacing is “*more likely related to under-resolved near-surface gradients and turbulent mixing at the boundary-layer top, to the [sub-grid scale] model formulation, and/or to numerical issues, and not to deficiencies due to the use of improper surface boundary conditions*”.

In the unstable case, the “FV free” scheme (with $\delta_a = z_{1\2}^{\text{low-res}}$) seems much more robust than the other schemes in the first 200 meters (remind that the height of the SL is approximately 10m). Above this height the differences between the high resolution and the low resolution simulations are not clearly influenced by the SL treatment. Note also that, as in the stable case, the “FV2” scheme is also less consistent: enforcing the MO theory in the first cell increases the sensitivity of the solution to $\delta_a$ because the SL is then tightly coupled with the computational domain. Finally, the “FV free” scheme combines good consistency properties with a SL scheme coherent with the physical theory.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Relative difference of the wind speed between low- and high-resolution simulations for several SL strategies. left: stable stratification. Right: unstable stratification.</figcaption></figure>

References

1. Clement, S.: Numerical analysis for a combined space-time discretization of air-sea exchanges and their parameterizations. PhD thesis, Université Grenoble Alpes (2022), ⟨tel-04066324⟩

2. Jaegle, F., Cabrit, O., Mendez, S., Poinsot, T.: Implementation methods of wall functions in cell-vertex numerical solvers. Flow Turbul Combust, 85, 245-272 (2010). doi:10.1007/s10494-010-9276-1

3. McWilliams, J. C., Huckle, E., and Shchepetkin, A. F.: Buoyancy Effects in a Stratified Ekman Layer. J. Phys. Oceanogr., 39, 2581-2599 (2009). doi:10.1175/2009JPO4130.1

4. Larsson, J., Kawai, S., Bodart, J., Bermejo-Moreno, I.: Large eddy simulation with modeled wall-stress: Recent progress and future directions. Mech. Eng. Rev., 3 (2016). doi:10.1299/mer.15-00418

5. Nishizawa, S., Kitamura, Y.: A Surface Flux Scheme Based on the Monin-Obukhov Similarity for Finite Volume Models. J. Adv. Model. Earth Syst., 12, 3159-3175 (2018). doi:10.1029/2018MS001534

6. Basu, S., Lacser, A.: A cautionary note on the use of Monin-Obukhov similarity theory in very high-resolution Large-Eddy Simulations. Bound.-Layer Meteorol., 163, pp. 351-355 (2017). doi:10.1007/s10546-016-0225-y

7. Piller, M. and Stalio, E.: Finite-volume compact schemes on staggered grids. Journal of Computational Physics, 197(1):299-340 (2004). doi:10.1016/j.jcp.2003.10.037

8. Clement, S.: Code for PhD thesis. Zenodo (2022). doi:10.5281/zenodo.7092357

9. Maronga, B., Knigge, C., Raasch, S.: An improved surface boundary condition for large-eddy simulations based on Monin-Obukhov similarity theory: evaluation and consequences for grid convergence in neutral and stable conditions. Bound.-Layer Meteorol., 174(2), 297-325 (2020). doi:10.1007/s10546-019-00485-w