# Designing a Fully-Tunable and Versatile TKE-1 Turbulence Parameterization for the Simulation of Stable Boundary Layers

Étienne Vignon, Khadija Arjdal, F Cheruy, M Coulon-Decorzens, C Dehondt, T Dubos, Sébastien Fromang, F. Hourdin, L Lange, L Raillard, et al.

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

## To cite this version:

Étienne Vignon, Khadija Arjdal, F Cheruy, M Coulon-Decorzens, C Dehondt, et al. "Designing a Fully-Tunable and Versatile TKE-1 Turbulence Parameterization for the Simulation of Stable Boundary Layers." Journal of Advances in Modeling Earth Systems, 2024, 10.1029/2024MS004400. hal-04737007

HAL Id: hal-04737007

https://hal.science/hal-04737007v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
RESEARCH ARTICLE

10.1029/2024MS004400

Key Points:

* A simple TKE-1 turbulent diffusion scheme is developed in a semi-heuristic way for applications in models of the Earth and Mars atmospheres

* All adjustable parameters are clearly identified and the number of parameters is minimized to thoroughly assess the parametric sensitivity

* Once tuned on GEWEX Atmospheric Boundary Layer Study 1 1D simulations, the scheme is able to capture the Antarctic and Martian stable boundary layers in 3D simulations

Correspondence to:

É. Vignon,
etienne.vignon@lmd.ipsl.fr

Citation:

Vignon, É., Arjdal, K., Cheruy, F., Coulon-Decorzens, M., Dehondt, C., Dubos, T., et al. (2024). Designing a fully-tunable and versatile TKE-1 turbulence parameterization for the simulation of stable boundary layers. *Journal of Advances in Modeling Earth Systems, 16*, e2024MS004400. https://doi.org/10.1029/2024MS004400

Received 15 APR 2024
Accepted 8 SEP 2024

© 2024 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Designing a Fully-Tunable and Versatile TKE-1 Turbulence Parameterization for the Simulation of Stable Boundary Layers

É. Vignon¹ ID, K. Arjdal¹,² ID, F. Cheruy¹ ID, M. Coulon-Decorzens¹, C. Dehondt³, T. Dubos¹ ID, S. Fromang³ ID, F. Hourdin¹ ID, L. Lange¹ ID, L. Raillard¹ ID, G. Rivière¹ ID, R. Roehrig⁴ ID, A. Sima¹ ID, A. Spiga¹ ID, and P. Tiengou¹,⁵

¹Laboratoire de Météorologie Dynamique-IPSL, Sorbonne Université/CNRS/ Ecole Normale Supérieure-PSL Université/ Ecole Polytechnique-Institut Polytechnique de Paris, Paris, France, ²International Water Research Institute (IWRI), CSAES — Mohammed VI Polytechnic University, Benguerir, Morocco, ³Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France, ⁴CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France, ⁵Milieux environnementaux, tranferts et interaction dans les hydrosystemes et les sols (Metis)/Sorbonne Université/IPSL/CNRS/ EPHE, Paris, France

**Abstract** This study presents the development of a so-called Turbulent Kinetic Energy (TKE)-1, or TKE-1, parameterization of the diffusion coefficients for the representation of turbulent diffusion in neutral and stable conditions in large-scale atmospheric models. The parameterization has been carefully designed to be completely tunable in the sense that all adjustable parameters have been clearly identified and the number of parameters has been minimized as much as possible to help the calibration and to thoroughly assess the parametric sensitivity. We choose a mixing length formulation that depends on both static stability and wind shear to cover the different regimes of stable boundary layers. We follow a heuristic approach for expressing the stability functions and turbulent Prandlt number in order to guarantee the versatility of the scheme and its applicability for planetary atmospheres composed of an ideal and perfect gas such as that of Earth and Mars. Particular attention has been paid to the numerical stability and convergence of the TKE equation at large time steps, an essential prerequisite for capturing stable boundary layers in General Circulation Models (GCMs). Tests, parametric sensitivity assessments and preliminary tuning are performed on single-column idealized simulations of the weakly stable boundary layer. The robustness and versatility of the scheme are assessed through its implementation in the Laboratoire de Météorologie Dynamique Zoom GCM and the Mars Planetary Climate Model and by running simulations of the Antarctic and Martian nocturnal boundary layers.

**Plain Language Summary** In planetary atmospheres, turbulent motions actively contribute to the mixing of quantities such as heat, momentum, and chemical species. Such motions are not resolved in coarse-grid atmospheric models and have to be parameterized. The parameterization of turbulent mixing should ideally be based on physical laws and sufficiently sophisticated to realistically represent the full spectrum of motions over the full range of stability encountered in the atmospheres. However, it also necessarily contains a number of closure parameters not always well identified and whose values are determined empirically—thereby questioning the universality of the parameterization and its potential application over the full globe or even to other planets—or adjusted to guarantee the numerical stability of the model. This study presents the design of a turbulent mixing parameterization that can be fully calibrated and applied in planetary atmospheres such as that of Mars. We then calibrate the parameterization on an idealized simulation set-up and test its robustness and performance by running simulations of the Antarctic and Martian atmospheres.

1. Introduction

Turbulence efficiently transports momentum, energy, moisture, and matter in the atmosphere, particularly in the planetary boundary layer where it controls sensible and latent heat fluxes as well as the transfer of momentum between the air and the ground surface. It thereby directly affects the diurnal cycle of the near-surface atmospheric quantities and also impacts on the lifetime and structure of synoptic-scale dynamical systems. Turbulent transport is therefore an essential component of the physics of climate models, numerical weather prediction models and more generally of General Circulation Models (GCMs) of planetary atmospheres. As turbulent eddies manifest on scales ranging from a few millimeters to a few tens of kilometers in deep convective systems, modelers develop
conceptually separated subgrid parameterizations targeting different types—or different scales—of transport processes. Non-local turbulent transport resulting from large and organized convective cells, being deep or shallow, is often treated with so-called mass flux schemes (e.g., Emanuel, 1991; Golaz et al., 2002; Hourdin et al., 2002; Tiedtke, 1989). Local turbulent mixing which results from eddies whose typical size is smaller or similar to the thickness of models' vertical layers—namely a few tens of meters—is often parameterized with a local K-gradient diffusion scheme. In those schemes, the turbulent flux is parameterized with a Fick's law type down-gradient diffusion formulation that relies on the introduction of a turbulent diffusion coefficients. Such schemes are particularly critical to simulate the stable and neutral atmospheric boundary layers (Cuxart et al., 2006; Delage, 1997; Sandu et al., 2013), the land-atmosphere coupling as well as the thermal inversion at the top of convective boundary layers.

Several K-gradient diffusion parameterizations have been developed since the pioneering work of Louis (1979) and have been the subject of a substantial body of literature in atmospheric sciences. Among them, the moderate-complexity 1.5 order schemes, or TKE-l schemes, express the diffusion coefficients as functions of a diagnostic vertical turbulent length-scale, or mixing length, and of a prognostic estimation of the Turbulent Kinetic Energy (TKE) (Mellor & Yamada, 1982; Yamada, 1983).

The closure of the TKE evolution equation and the empirical and/or heuristic formulation of the mixing-length necessarily introduce free parameters in the parameterization, and therefore a certain degree of empiricism in the expression of the diffusion coefficients (Li et al., 2016). Indeed, such parameters do not have, by essence, fixed and universal values. Some of them—and the associated variability range thereof—are determined empirically using field observations, laboratory data, Large Eddy Simulations (LES) or Direct Numerical Simulations (DNS) while others are arbitrarily set. In practice, in climate and numerical weather prediction models, the value of some coefficients is often retuned to match large-scale or meteorological targets. For instance as all mixing processes that occur at smaller scales than typical GCM mesh scale are not parameterized—such as small scale internal waves or submeso-scale motions—the mixing in stable conditions is often artificially enhanced to prevent unrealistic runaway surface cooling due to surface-atmosphere mechanical decoupling and to maintain sufficient surface drag and Ekman pumping in extratropical cyclones (Holtslag et al., 2013; Sandu et al., 2013). Such empiricism and Earth-oriented tuning can somewhat question the applicability of these turbulent mixing parameterizations in planetary GCMs, even in GCMs of Mars (e.g., Colaïtis et al., 2013; Forget et al., 1999) where the planetary boundary layer shares similarities with that on Earth (Spiga et al., 2010).

In addition, arbitrary parameter calibration—sometimes beyond reasonable ranges—is often required to improve the numerical convergence and stability of the parameterization once it is implemented in models with typical physics time steps of a few minutes to a few tens of minutes. Indeed, the numerical implementation of a K-gradient turbulence scheme is prone to spurious oscillations called “fibrillations” (Girard & Delage, 1990; Kalnay & Kanamitsu, 1988). Such fibrillations are due to (a) the coupling between momentum and potential temperature via the turbulent diffusion coefficients and (b) the discretization of the vertical diffusion in which the nonlinear exchange coefficient is often treated explicitly in time. Even though the TKE budget is often close to a local equilibrium (Lenderink & Holtslag, 2004), the prognostic prediction of the TKE generally makes TKE-l schemes less sensitive to the time discretization and less prone to fibrillation than traditional first-order schemes (Bazile et al., 2011) in which the diffusion coefficients are explicit and diagnostic functions of the mean static stability and wind shear (Delage, 1997; Louis, 1979; Louis et al., 1982). This is mostly explained by the fact that the prognostic TKE plays a role of “reservoir” that damps the sometimes abrupt evolution of the diffusion coefficients with time (Mašek et al., 2022). However, even TKE-based schemes can be affected by numerical instabilities which can be related to the numerical treatment of the TKE equation itself (Deleer-snijder, 1992; Vignon et al., 2018) or to the coupling with other prognostic quantities such as the turbulent potential energy (Mašek et al., 2022). The numerical treatment of the TKE equation and more generally of the turbulent diffusion thereby come out as a forefront issue in atmospheric modeling. Hence, one has to find a good trade-off between the complexity and sophistication of a turbulent mixing scheme and its practical implementation in large scale atmospheric models avoiding as much as possible unrealistic parameter calibration to guarantee numerical stability and fair model performances.

The sensitivity of the stable boundary layer representation to turbulent diffusion calibration in a large scale atmospheric model was assessed in a game-changing study by Audouin et al. (2021) using a semi-automatic tuning tool based on uncertainty quantification (Couvreux et al., 2021; Hourdin et al., 2021). The authors identified a few
key tuning parameters—and their acceptable ranges of values—in the TKE-1 turbulent diffusion scheme of the Arpege-Climat model and assessed to what extent biases in the simulation of the extremely stable Antarctic boundary layer are explained by structural parameterization deficiencies or tuning choices. However, the boundary layer and surface layer schemes of Arpege-Climat contain a large number of tuning parameters, sometimes subtly interdependent, and considering all of them in a tuning exercise may be confusing, thereby challenging.

The present study aims to design a new and simple TKE-1 turbulent diffusion scheme for large scale atmospheric models

1. that is sufficiently robust and versatile to be applicable on both Earth and Mars, and potentially on other planetary atmospheres and;

2. that is built to be completely tuned in the sense that all adjustable parameters are clearly identified and their number minimized to help the calibration—or parameter adjustment—and assess the parametric sensitivity.

The scheme will be referred to as the ATKE scheme—for Adjustable TKE-1 scheme—in the paper.

We follow a simple heuristic approach—as in Lenderink and Holtslag (2004) and He et al. (2019)—for expressing the stability functions and turbulent Prandlt number to guarantee the versatility of the scheme and its potential applicability for planetary atmospheres composed of an ideal and perfect gas. A particular attention is paid to the numerical treatment of the TKE prognostic equation to ensure the numerical stability even in conditions of strong wind shear or strong stratification. It is worth emphasizing that the “local” nature of the scheme makes it mostly appropriate for neutral and stably stratified conditions, hence the particular focus on stable boundary layers in the paper. The scheme is tested and tuned—using the same Uncertainty Quantification approach as in Audouin et al. (2021) and Hourdin et al. (2021)—on idealized single column simulations of the stable boundary layer. The parameterization is then implemented and tested in the Earth Laboratoire de Météorologie Dynamique Zoom (LMDZ) GCM (Cheruy et al., 2020; Hourdin et al., 2020) and the Mars Planetary Climate model (Forget et al., 1999) to verify its robustness and assess its performances in simulating the stable Antarctic and Martian nocturnal boundary layers.

## 2. Parameterization Development

This section presents the derivation of the ATKE scheme, starting briefly and purposely with some generalities to clearly set the parameterization in the framework of turbulent diffusion in GCMs of planetary atmospheres.

### 2.1. General Framework

The conservation law for an extensive quantity c—being for example, the potential temperature, wind components or concentration in chemical species—in a compressible atmosphere reads:

$$ \frac{\partial \rho c}{\partial t} + \vec{\nabla}(\rho \tilde{u} c) = P_c \quad (1) $$

With, in Cartesian coordinates $(x, y, z)$, $\tilde{u} = \tilde{u} + \tilde{v} + \tilde{w}$ the wind vector, $\rho$ the air density and $P_c$ the net source/loss term. We note the statistical (ensemble) average with an overline and introduce the air weighting average operator $\sim$ such that:

$$ \tilde{c} \equiv \frac{\overline{\rho c}}{\overline{\rho}} \quad (2) $$

Note that $\tilde{c}$ is an extensive physical quantity per mass unit. We decompose c into a mean state and a fluctuation such that $c = \tilde{c} + c'$. We then apply the statistical average operator (overline) on Equation 1 that now reads:

$$ \underbrace{\frac{\partial \overline{\rho \tilde{c}}}{\partial t} + \overline{\nabla}(\rho \tilde{c} \tilde{u}}}_{(1)} = - \underbrace{\overline{\nabla}(\overline{\rho \tilde{u}' c'}}}_{(2)} + \overline{P_c} \quad (3) $$
In large-scale atmospheric models the scale separation is imposed by the size of the grid cells which determines the resolved and unresolved components. In this framework, the term (1) in Equation 3 is handled by the dynamical core while the term (2) is the essence of the physical subgrid parameterizations. Further assuming that the horizontal fluxes can be ignored compared to vertical ones given the typical horizontal ($\Theta(10-100 \text{ km})$) and vertical ($\Theta(10-100 \text{ m})$) resolutions of atmospheric models, it follows that $\overline{\nabla}(\overline{\rho u'}c') \approx \partial_z(\overline{\rho w'}c')$. A local turbulent mixing parameterization aims at calculating a tendency on the mean state variable $\tilde{c}$ due to the vertical turbulent diffusion as follows:

$$ \left. \frac{\partial \tilde{c}}{\partial t} \right|_{\text{diffusion}} = - \frac{1}{\bar{\rho}} \frac{\partial \overline{\rho w' c'}}{\partial z} \quad (4) $$

For better readability and conciseness, we leave the ~ notation out for mean state quantities and note $\rho = \bar{\rho}$ in the following.

For local and mostly shear driven turbulent eddies, the mixing of any conservative quantity during turbulent mixing—such as the common Betts (1973)' variables—can be represented as a diffusive process (e.g., Louis, 1979). Turbulent fluxes can then be expressed with a down-gradient form: $\overline{\rho w' c'} = -\rho K_c \partial_z c$, $K_c$ being a diffusion coefficient. Equation 4 hence reads:

$$ \left. \frac{\partial c}{\partial t} \right|_{\text{diffusion}} = \frac{1}{\bar{\rho}} \frac{\partial}{\partial z} \left( \rho K_c \frac{\partial c}{\partial z} \right) \quad (5) $$

Once the $K_c$ coefficient has been calculated at vertical model layer interfaces, such an equation can be numerically solved with an implicit approach through the inversion of a tri-diagonal matrix.

We now focus on the closure of the $K_c$ coefficient which is the main scope of the present study. We follow here an approach historically proposed by Mellor and Yamada (1974) and Yamada and Mellor (1975) that is, a 1.5 order closure or TKE-l scheme. In this framework, $K_c$ coefficient is expressed as the product of a vertical turbulent length scale or mixing length $l$ with a turbulent vertical velocity scale taken proportional to the square root of the TKE $e = \frac{1}{2}(u'^2 + v'^2 + w'^2)$. The latter is multiplied by a stability function $S_c$ that accounts for the fact that the turbulence anisotropy—thus the contribution of TKE to vertical turbulent mixing—varies with the local stability of the atmosphere characterized by the gradient Richardson number $Ri$. The diffusion coefficient $K_c$ is then expressed as (Yamada, 1983; Zilitinkevich et al., 2007):

$$ K_c = l S_c (Ri) \sqrt{e} \quad (6) $$

In the following sections, we describe the estimation of the three different terms of $K_c$, namely $e$, $S_c$, and $l$. As we want our turbulent scheme to be applicable on Earth and Mars (and potentially other planetary environments), we have to ensure that their expressions are as planet-independent as possible.

## 2.2. TKE Prognostic Equation

### 2.2.1. Parameterization of the Source and Loss Terms

Assuming the horizontal homogeneity of the subgrid-scale statistics, the TKE obeys the following evolution equation (Stull, 1990):

$$ \frac{\partial e}{\partial t} = \underbrace{-\overline{u'w' \frac{\partial u}{\partial z} - v'w' \frac{\partial v}{\partial z}}}_{W} + \underbrace{\overline{b'w' \frac{1}{\rho} \frac{\partial}{\partial z}(\overline{\rho w' e} + \overline{w' p})}}_{B} - \underbrace{\overline{-\epsilon}}_{T} \quad (7) $$

$W$ is the wind shear production term that can be expressed with the down-gradient expression of fluxes with a diffusion coefficient for momentum hereafter denoted as $K_m$:
$$-\overline{u'w'} \frac{\partial u}{\partial z} - \overline{v'w'} \frac{\partial v}{\partial z} = K_m S^2 = lS_m \sqrt{e} S^2 \quad (8)$$

with $S^2 = (\partial_z u)^2 + (\partial_z v)^2$ the wind shear and $S_m$ the stability function for momentum. $\mathcal{B}$ is the buoyancy b production/consumption term. For dry air under the ideal gas assumption, one can write:

$$\overline{b'w'} = \frac{-g\partial\rho}{\rho\partial\theta}|_p \overline{w'\theta'} = \frac{g}{\theta}\overline{w'\theta'} = -K_h \frac{g}{\theta}\frac{\partial\theta}{\partial z} = -K_h N^2 = -lS_h\sqrt{e}N^2 \quad (9)$$

where $g$ is the gravity acceleration of the planet, $\theta$ the potential temperature, $N$ the Brünt-Väisälä pulsation, $K_h$ the diffusion coefficient for heat and $S_h$ the stability function for heat. In the case of an atmosphere containing water vapor or chemical species $\xi$, buoyancy reads $\overline{b'w'} = \frac{-g}{\rho}(\frac{\partial\rho}{\partial\theta}|_{p,\xi}\overline{w'\theta'} + \frac{\partial\rho}{\partial\xi}|_{p,\theta}\overline{w'\xi'})$. For water vapor—in absence of phase change—or for non-reactive chemical species, one can define a virtual temperature $T_v$ (and a subsequent virtual potential temperature $\theta_v$) corresponding to the temperature that dry air would have if its pressure and density were equal to those of a given sample of the mixture of gas. In this case:

$$\overline{b'w'} \approx \frac{g}{\theta_v}\overline{w'\theta_v} = -\frac{g}{\theta_v}K_h\frac{\partial\theta_v}{\partial z} \quad (10)$$

It is worth noting here that the expression of the buoyancy term (or Brünt-Väisälä pulsation) is gravity-dependent thus planet-dependent. For simplicity and consistency with previous literature on turbulent mixing schemes, we keep the formalism with explicit gravity in the following. However, a more universal derivation of the scheme can be achieved with a gravity-invariant formulation of the TKE and turbulent diffusion equations. Such a formulation is proposed in Appendix A.

$\mathcal{D}$ is the viscous TKE dissipation term that can be expressed following Kolmogorov (1941):

$$\epsilon = \frac{e^{3/2}}{l_c} \quad (11)$$

with $l_c$ the dissipation length-scale characterizing the size of the most dissipative and energy-containing eddies. Following for instance Yamada (1983) and Bougeault and Lacarrère (1989), we assume that $l_c$ scales with $l$ such that $l_c = c_e l$, $c_e$ being a scalar. Its value roughly ranges between 1.2 and 10.0 (Audouin et al., 2021; He et al., 2019; Yamada, 1983) since dissipation length scale—characterizing the dissipation of turbulence as a whole—might be larger than vertical mixing length in stable conditions due to the fact that kinetic energy can dissipate through wavy motion with little transfer to the smaller turbulent scales (Cuxart et al., 2006).

The vertical turbulent flux of TKE and the pressure term gathered in $\mathcal{T}$ redistribute TKE through the depth of the atmospheric column. Hence, those two terms are commonly grouped together and expressed as a TKE turbulent diffusion term:

$$-\frac{1}{\rho}\frac{\partial}{\partial z}(\overline{\rho w'}e + \overline{w'p'}) = \frac{1}{\rho}\frac{\partial}{\partial z}(\rho K_e \frac{\partial e}{\partial z}) \quad (12)$$

$K_e$ being taken proportional to $K_m$ (Bougeault & Lacarrère, 1989; Lenderink & Holtslag, 2004; Yamada, 1983): $K_e = c_e K_m$. $c_e$ is a constant whose value is generally around 1–2 and that we will arbitrarily allow to vary between 1 and 5 (Baas et al., 2018; Bougeault & Lacarrère, 1989; Lenderink & Holtslag, 2004). The lower boundary condition of $e$ that is, the surface value of the TKE $e_s$, is estimated by assuming stationary near-neutral conditions in the surface layer. On such a condition (Baas et al., 2018; Lenderink & Holtslag, 2004):

$$e_s = c_s u_*^2 \quad (13)$$
with $c_s$ a constant and $u_*$ the surface friction velocity calculated from the surface drag coefficient for momentum and the wind speed at the first model level. A proper scaling of the TKE-l parameterization with the Monin-Obukhov similarity in the surface layer requires (He et al., 2019):

$$c_s = c_e^{2/3} \qquad (14)$$

### 2.2.2. Numerical Treatment

Once the different TKE source and loss terms have been expressed, Equation 7 has to be integrated in time. The numerical treatment of Equation 7 is critical as the solution must be stable and converge at typical physical time steps used in atmospheric GCMs namely, of the order of ≈15 min. Several methods have been proposed in the literature, particularly regarding the treatment of the dissipation term with different degrees of implicitness (Bazile et al., 2011).

Here, we propose a two-step resolution method which allows for an exact treatment of the dissipation term under some assumptions, while the transport term is calculated separately.

*Step 1* We calculate the TKE tendency due to the shear, buoyancy and dissipation terms. Noting $q = \sqrt{2e}$, one can rewrite Equation 7 with no transport term as:

$$\frac{\partial q}{\partial t} = \frac{IS_m}{\sqrt{2}} S^2 \left(1 - \frac{Ri}{Pr}\right) - \frac{q^2}{2^{3/2}c_e l} \qquad (15)$$

with $Pr = \frac{K_m}{K_h} = \frac{S_m}{S_h}$ the turbulent Prandtl number. We then solve this equation through an implicit treatment of $q$ assuming that the mean temperature and wind field does not vary much during the time step $\delta t$ and thus keeping the explicit value—that is the value at the beginning of the time step—of $Ri$, $S_m$, $Pr$, and $l$. Equation 16 then reads:

$$\frac{q_{t+\delta t} - q_t}{\delta t} = \frac{IS_m}{\sqrt{2}} S^2 \left(1 - \frac{Ri}{Pr}\right) - \frac{q_{t+\delta t}^2}{2^{3/2}c_e l} \qquad (16)$$

then can be rewritten in a second-order polynomial form after some rearrangement:

$$q_{t+\delta t}^2 + A_t q_{t+\delta t} + B_t = 0 \qquad (17)$$

$$\text{with } A_t = \frac{c_e l^{2/3}}{\delta t} \text{ and } B_t = -\left(\frac{q_t c_e l^{2/3}}{\delta t} + 2l^2 c_e S_m S^2 \left(1 - \frac{Ri}{Pr}\right)\right)$$

One can show that given the choice we will make for the formulation of the turbulent Prandtl number in the next section, $Ri/Pr$ namely the flux Richardson number, is by construction always <1. This in fact reflects a condition imposed by steady-state TKE budget equation for which the wind shear production term and the buoyancy term cannot exceed unity to maintain a non-zero TKE dissipation thus a non-zero turbulence (e.g., Zilitinkevich et al., 2008).

The discriminant $\Delta = A_t^2 - 4B_t$ of Equation 17 is thus always >0 and Equation 17 always admits a positive solution for $q$ thus $e$ that reads:

$$e = \frac{(-A_t + \sqrt{\Delta})^2}{8} \qquad (18)$$

*Step 2* The TKE variation due to the transport term $\mathcal{T}$ is then calculated and added to the value found in step 1. The calculation of this term consists in resolving the following equation:
<figure><img src="image_4.png" /><figcaption>Figure 1. $S_{m,h}$ (panel a) and Pr (panel b) as functions of the Richardson number $Ri$ following Equations 20 and 23. Envelopes show the range of variation when adjustable parameters evolve in their range of acceptable values (Table 1). Solid lines show the curves for the following arbitrary set of parameters' values: $c_e = 5.9$, $Pr_n = 0.8$, $\alpha_{Pr} = 4.5$, $r_\infty = 2$, $Pr_\infty = 0.4$, $S_{min} = 0.05$ and $Ri_c = 0.2$.</figcaption></figure>

$$ \frac{\partial e}{\partial t} = \frac{1}{\rho} \frac{\partial}{\partial z} \left( \rho K_e \frac{\partial e}{\partial z} \right) \qquad (19) $$

with an a priori knowledge of $K_e$—namely an explicit value of $K_e$ calculated with the $e$ value from Step 1—Equation 19 is a typical diffusion equation that is solved implicitly in time through a tri-diagonal matrix inversion (Dufresne & Ghattas, 2009).

We will see hereafter that this proposed two-step method performs well in terms of numerical stability and convergence. A brief comparison (not shown) with the non-split TKE equation resolution method used for Météo France models (Marquet, 2008)—with a semi-implicit treatment of the dissipation term—exhibits quite similar numerical performances. Our method appears slightly more stable but also slightly more diffusive at large time steps, which is likely due to the full-implicit treatment of the TKE source and loss terms.

## 2.3. Heuristic Expressions for the Stability Functions and Turbulent Prandtl Number

We now have to derive a heuristic expression for the stability function $S_m$ of the gradient Richardson number $Ri = N^2/S^2$ to be used in the formulation of the diffusion coefficient for momentum. On one hand, $S_m$ should increase when an atmospheric layer locally becomes more unstable and thus with decreasing negative $Ri$. On the other hand, we want to prevent $S_m$ from reaching infinite value when $Ri \to -\infty$ to avoid risk of numerical instabilities when $K_m \to \infty$ (Lenderink & Holtslag, 2000). It is worth recalling here that in unstable conditions, turbulent transport becomes non-local and another type of parameterization such as a mass-flux scheme should come in support of the K-diffusion. In stable conditions as turbulent mixing intensity decreases with increasing stability, we assume a simple linear decrease with $Ri$ down to a minimum value attained when the Richardson number equals a critical value (Mellor & Yamada, 1974).

Following Lenderink and Holtslag (2004), we propose the following expression for $S_m$ plotted in Figure 1a:

$$ S_m(Ri) = \begin{cases} c_n + \frac{2}{\pi}(c_\infty - c_n) \arctan\left(\frac{-Ri}{Ri_0}\right) & \text{if } Ri < 0 \\ \max\left(c_n\left(1 - \frac{Ri}{Ri_c}\right), S_{min}\right) & \text{if } Ri \ge 0 \end{cases} \qquad (20) $$

$c_n$ is the value of $S_m$ at $Ri = 0$ and $c_\infty$ is the $S_m$ value in the convective limit. $r_\infty = c_\infty/c_n$ is comprised between 1.2 and 5 (Lenderink & Holtslag, 2004; Mellor & Yamada, 1982). $Ri_c$ is a critical Richardson number whose inverse value controls the slope of $S_m$ in stable conditions. Previous literature suggests $Ri_c$ values comprised between 0.19
and 0.25 (He et al., 2019; Mellor & Yamada, 1974, 1982). As the turbulence vertical anisotropy, that is $\frac{1}{2}w^2/e$, does not reach 0 in very stable conditions (Li et al., 2016; Zilitinkevich et al., 2007), $S_m$ must be lower-bounded by a value $S_{min}$ which is roughly around 0.05 and that we will make vary between 0.025 and 0.1.

The continuity in slope for $Ri = 0$ further gives:

$$Ri_0 = \frac{2}{\pi}(c_\infty - c_n) \frac{Ri_c}{c_n} \qquad (21)$$

Furthermore, the so-called local-scaling similarity theory in stable boundary layers (Derbyshire, 1990; Nieuwtsadt, 1984; van de Wiel et al., 2010) implies that in stationary conditions, turbulent fluxes and vertical gradient wind speed must scale such that $\frac{K_m}{l^2_S}$ converges toward 1 in the neutral limit. This conditions leads to a direct relationship between $c_n$ and the coefficient $c_c$ (Baas et al., 2018; He et al., 2019), the latter being the ratio between the mixing length $l$ and the TKE dissipation length scale (Section 2.2.1):

$$c_n = c_c^{-1/3} \qquad (22)$$

The stability function for the heat flux $S_h$ is estimated through a parametrization of the turbulent Prandtl number Pr. Under unstable conditions, the dominant coherent structures such as rising plumes and thermals have vertical velocity anomalies which generally better correlate with buoyancy and temperature anomalies than momentum anomalies in average. Therefore, one expects Pr to decrease with increasing instability (Li, 2019). In stably stratified conditions, buoyancy is expected to suppress the transport of heat but the existence of gravity waves can maintain some transport of momentum inducing an increase in Pr with increasing stability. Collection of field experiments, laboratory data, LES and DNS results show a consistent increase in Pr with $Ri$ with a asymptotical linear behavior at strong stability (Li, 2019; Zilitinkevich et al., 2008). We therefore propose the following expression of Pr that is plotted in Figure 1b:

$$Pr(Ri) = \begin{cases} \frac{Pr_n - \frac{2}{\pi}(Pr_\infty - Pr_n)}{\arctan(\frac{-Ri}{Ri_1})} & \text{if } Ri < 0 \\ Pr_n e^{\frac{1-\alpha Pr_n}{Pr_n}Ri} + \alpha_{Pr}Ri & \text{if } Ri \ge 0 \end{cases} \qquad (23)$$

The formulation in stable conditions is inspired from Venayagamoorthy and Stretch (2010) and it shows fair agreement with experimental data (Li, 2019). $\alpha_{Pr}$ is the slope of the asymptotical linear trend at high stability and its value ranges from 3 to 5 (Grisogono, 2010). $Pr_n$ is the neutral value of Prandtl number which from extensive laboratory and field experiments as well as theoretical works range from 0.7 to 1 (Grisogono, 2010; Li, 2019). The continuity in slope at $Ri = 0$ gives

$$Ri_1 = \frac{2}{\pi}(Pr_\infty - Pr_n) \qquad (24)$$

$Pr_\infty$ is the value of Pr in the convective limit and its value roughly ranges between 0.3 and 0.5 (Li, 2019).

## 2.4. Vertical Turbulent Mixing Length Formulation

In near-neutral conditions, we choose a turbulent vertical length-scale formulation $l_n$ similar to Blackadar (1962) in which the displacement of eddies is limited by the distance to the ground in the neutral limit:

$$l_n = \frac{\kappa z l_\infty}{\kappa z + l_\infty} \qquad (25)$$

where $\kappa$ is the Von Kármán constant. $l_\infty$ is the mixing-length far above the ground whose value in near-neutral conditions is generally estimated between 15 and 75 m (Lenderink & Holtslag, 2004; Sun, 2011). In stable conditions, the vertical displacement of eddies—whose size is roughly above the so-called Ozmidov scale—is limited by the stratification of the flow (e.g., van de Wiel et al. (2008)). André et al. (1978) and
Deardoff (1980) introduced a widely used buoyancy length-scale which depends on the flow stratification characterized by Brunt-Väisälä pulsation N. The mixing length in stable conditions $l_s$ then read:

$$l_s = c_l \frac{\sqrt{e}}{N} \qquad (26)$$

$c_l$ being a scalar whose value varies between 0.1 and 2 (Baas et al., 2018; Deardoff, 1980; Grisogono & Belušić, 2008; Nieuwtsadt, 1984).

More recent studies introduced wind-shear dependent formulation of $l_s$ to account for the deformation of eddies— whose size is above a so-called Corrsin scale—by vertical wind shear (e.g., Grisogono, 2010; Grisogono & Belušić, 2008; Rodier et al., 2017). Grisogono and Belušić (2008) proposed a mixing-length formulation including both the effect of stratification and vertical wind shear $S^2$ that reads:

$$l_s = c_l \frac{\sqrt{e}}{2\sqrt{S^2}(1 + \sqrt{Ri}/2)} \qquad (27)$$

The final mixing-length $l$, being either ground-limited or stratification-limited is the minimum between $l_n$ and $l_s$. In the model implementation, we choose a commonly-used continuous inverse interpolation formulation:

$$l = \left( \frac{1}{l_n^\delta} + \frac{1}{l_s^\delta} \right)^{-1/\delta} \qquad (28)$$

$\delta$ is set to 1 by default to guarantee a smooth transition between the neutral and the stratification length-scales that is, a smooth transition between the so-called “wall” effect of the ground and the buoyancy and wind shear effects. The two expressions of $l_s$ can be used independently in the parameterization but unless otherwise stated, the results presented in the rest of the paper have been obtained with formulation dependent on both stratification and wind shear (Equation 27). As $l_s$ depends on the TKE, $l$ is calculated with an explicit value of the TKE that is, the value at the beginning of the time-step. In practice, $l$ is also lower bounded by a value $l_{min} = 1$ cm to prevent it from reaching value below the Kolmogorov length scale in planetary atmospheric motions (Chen et al., 2016). Such a minimum value is also necessary to allow the TKE to have a weak but non-null steady-state limit value at very high stability that is, at very high $Ri$ numbers.

## 2.5. Surface Layer Scheme Matching

In line with England and McNider (1995) and Redelsperger et al. (2001), we propose here a derivation of the surface drag coefficients that is consistent with the ATKE parameterization by matching the ATKE scheme—in stationary conditions—with the surface similarity theory.

Neglecting the vertical diffusion term of TKE that is, $\mathcal{T}$, Equation 7 in stationary conditions ($\partial_t e = 0$) can read:

$$0 = K_m S^2 \left(1 - \frac{Ri}{P_r}\right) - \frac{e^{3/2}}{c_e l} = K_m S^2 \left(1 - \frac{Ri}{P_r}\right) - \frac{K_m^3}{c_e l^4 S_m^3} \qquad (29)$$

It can then be rearranged to give a first-order turbulent closure like expressions of the eddy diffusion coefficients for momentum and heat (Cuxart et al., 2006):

$$K_m = l^2 \sqrt{S^2} F_m(Ri) \qquad (30)$$

$$K_h = K_m / P_r = l^2 \sqrt{S^2} F_h(Ri) \qquad (31)$$

where
<table><caption>**Table 1**
Name, Definition, and Range of Acceptable Values for the Adjustable Parameters</caption><thead><tr><th>Name</th><th>Definition</th><th>Range</th></tr></thead><tbody><tr><td>c<sub>e</sub></td><td>Controls the value of the dissipation length scale</td><td>[1.2–10]</td></tr><tr><td>c<sub>e</sub></td><td>Controls the value of the diffusion coefficient of TKE</td><td>[1–5]</td></tr><tr><td>l<sub>∞</sub></td><td>Asymptotic mixing length far from the ground</td><td>[15–75]</td></tr><tr><td>c<sub>l</sub></td><td>Controls the value of the mixing length in stratified conditions</td><td>[0.1–2]</td></tr><tr><td>Ri<sub>c</sub></td><td>Critical Richardson number controlling the slope of S<sub>m</sub> in stable conditions</td><td>[0.19–0.25]</td></tr><tr><td>S<sub>min</sub></td><td>Minimum value of S<sub>m</sub> in very stable conditions</td><td>[0.025–0.1]</td></tr><tr><td>Pr<sub>n</sub></td><td>Neutral value of the Prandtl number</td><td>[0.7–1]</td></tr><tr><td>α<sub>Pr</sub></td><td>Linear slope of Pr with Ri in the very stable regime</td><td>[3–5]</td></tr><tr><td>r<sub>∞</sub></td><td>Ratio between c<sub>∞</sub> and c<sub>n</sub> controlling the convective limit of S<sub>m</sub></td><td>[1.2–5.0]</td></tr><tr><td>Pr<sub>∞</sub></td><td>Value of Pr in the convective limit</td><td>[0.3–0.5]</td></tr></tbody></table>

Note. Parameters are dimensionless exception $l_∞$ which is a length in m. Parameters in bold are those which affect the simulation of the neutral and stable boundary layer.

$$F_m(Ri) = S_m^{3/2} \sqrt{c_e} \left(1 - \frac{Ri}{Pr}\right)^{1/2} \qquad (32)$$

$$F_h(Ri) = S_m^{7/4} Pr^{-1} \sqrt{c_e} \left(1 - \frac{Ri}{Pr}\right)^{1/2} \qquad (33)$$

are first-order like stability functions. Near the ground in the surface layer, $l \approx \kappa z$ and England and McNider (1995) then show that $F_{m,h}$ functions are identical to the stability functions involved in the bulk expressions of the surface drag coefficients used to calculate surface fluxes of momentum and heat $C_{m,h}$ in models:

$$C_{m,h} = \frac{\kappa^2}{\log(z/z_{0m}) \log(z/z_{0m,h})} F_{m,h} \qquad (34)$$

with $z_{0m}$ and $z_{0h}$ the surface roughness lengths for momentum and heat respectively. Provided turbulence in the surface layer can be assumed to be close to a stationary state, using the same formulations for $S_m$ and $Pr$ in both the turbulent diffusion and surface layer schemes leads to a fully consistent formulation of turbulent fluxes from the surface layer up to the top of the boundary-layer.

## 2.6. Degrees of Freedom of the Scheme and Adjustable Parameters

Table 1 summarizes all the 10 adjustable parameters of the new parameterization and their ranges of acceptable values as previously introduced in the text. The eight first parameters in bold are those affecting the simulation of the neutral and stable boundary layers and taken into account in the tuning phase in the next section. It is worth mentioning that we also lower-bound the turbulent diffusion coefficients with the kinematic molecular viscosity and conductivity of the air, which are not tuning parameters per se but pressure and temperature dependent—thus planet dependent—quantities.

# 3. Implementation in General Circulation Models, Evaluation and Tuning

## 3.1. Implementation in the LMDZ GCM and Mars Planetary Climate Model

The ATKE parameterization has been implemented in the LMDZ Earth GCM (Cheruy et al., 2020; Hourdin et al., 2020), atmospheric component of the French IPSL Coupled-Model (Boucher et al., 2020) involved in the Coupled Model Intercomparison Project exercises. The turbulent-mixing parameterization of LMDZ has received a lot of attention in the past two decades, particularly regarding the convective boundary layer and the very stable boundary layer. It is a hybrid scheme in the sense that in convective boundary layers, turbulent fluxes are
expressed as a sum of a K-diffusion term—from the TKE-1 scheme of Yamada (1983) and revisited in Hourdin et al. (2002) and Vignon, Hourdin, et al. (2017)—and a non-local transport term by convective plumes parameterized with a mass-flux scheme from Hourdin et al. (2002), Rio et al. (2010), Hourdin et al. (2019). In this case, the turbulent flux of a quantity c reads:

$$ \overline{\rho w' c'} = -\rho K_c \frac{\partial c}{\partial z} + \rho \dot{\hat{\alpha}}(\hat{c} - c) \quad (35) $$

where $\hat{\alpha}$ is the fraction of the horizontal surface covered by ascending plumes, $\hat{w}$ the turbulent vertical velocity in the thermals, $\hat{c}$ the value of c in the thermals.

Despite the development efforts of the previous years, recent tests revealed that the latest version of the model—the CMIP6 version—still exhibits numerical instabilities in near-neutral boundary layers in presence of strong wind shear. In our implementation, the ATKE scheme exactly replaces the previous TKE-1 scheme which treats the turbulent diffusion in the whole atmosphere, with no specific triggering criterion. Similarly, in convective boundary layers the ATKE scheme works together with the LMDZ mass-flux scheme.

As a proof of concept, the ATKE scheme has also been implemented in the Mars Planetary Climate Model (Mars PCM, Forget et al. (1999)). This model also uses a hybrid scheme with a TKE-1 diffusion scheme inspired from Yamada (1983) and a dry parameterization of convective plumes (Colaïtis et al., 2013). Colaïtis et al. (2013) have pointed out that the default TKE-1 scheme of Hourdin et al. (2002) leads to numerical oscillations in strongly stratified Martian nighttime conditions. They addressed this issue by imposing a minimum mixing coefficient $K_{min}$ whose value depends on the boundary layer height following Holtslag and Boville (1993).

## 3.2. Parametric Sensitivity of the ATKE Scheme and Tuning

### 3.2.1. Initial Test on the GABLS1 Case and Parametric Sensitivity

The ATKE scheme is first tested on single column simulations using the 1D version of LMDZ. We use the 95-level vertical grid introduced in Hourdin et al. (2019) with a first model level at about 10 m above the ground and a first layer interface—where the eddy diffusivity coefficients are defined—at about 20 m above the ground. We run 1D simulations on the GEWEX Atmospheric Boundary Layer Study 1 (GABLS1) single column model intercomparison exercise. This exercise is often used to evaluate the performance of atmospheric turbulent diffusion schemes (e.g., Lopez-Gomez et al. (2020)) and consists in a no-radiation idealized 9 hr simulation of the development of a weakly stable boundary layer, with a constant zonal geostrophic wind of 8 m s⁻¹ and a constant surface cooling of −0.25 K h⁻¹ (Cuxart et al., 2006). The surface-layer Richardson number varies from 0 to ≈ 0.028 during the run. The fair convergence of 3D LES on this case—with the exact same initial and boundary conditions as those for single column models—make LES suitable references for GABLS1. Nonetheless, to sample the small variability between LES runs, we consider hereafter five reference LES which correspond to the MO-1m, MO-2m, UIB-2m, IMUK-1m, IMUK-2m simulations listed in Table 2 of Beare et al. (2006), the suffix referring to the vertical resolution.

Given the ranges of acceptable values associated with each of the n = 8 free parameters affecting the simulation of the stable boundary layer listed in Table 1, we need to run simulations with different sets of parameters to assess the parametric sensitivity of the scheme. For this purpose, we use the HighTune explorer statistical tool originally developed in the Uncertainty Quantification community and now applicable in atmospheric modeling (Couvreux et al., 2021). This tool allows to make a first perturbed physics ensemble experiment through an exploration of the initial n-dimension hypercube of parameters defined by the intervals given in Table 1 using a Latin Hyper Cube sampling method. Here 80 (10 times n) sets of parameters or free parameters' vectors are sampled. Unless otherwise stated, the simulations are run with a 15 min time step, that is, the typical value used for the LMDZ physics and that used for the ensemble of CMIP6 simulations.

Figure 2 shows the results of this a priori sensitivity analysis to free parameters' values for the vertical profiles of potential temperature, wind speed and TKE averaged over the eighth hour of the simulation. The yellow envelope displays the variability (minimum and maximum values) amongst the 80 simulations from this first so-called “wave” of simulations. Here, a wave is an ensemble of 80 simulations for which each simulation corresponds to a given vector of parameters sampled in the Latin Hypercube. Albeit encompassing the five reference LES
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Evolution of envelopes of the vertical profiles of potential temperature (panel a), wind speed (panel b) and Turbulent Kinetic Energy (TKE) (panel c) after 9 hr of GEWEX Atmospheric Boundary Layer Study 1 simulation. Yellow and orange envelopes correspond to waves 1 and 20 respectively that is, to the 1st and 20th set of 80 simulations during the tuning exercise. Blue curves show the five reference Large Eddy Simulations (LES). The red curve shows the “best” LMDZ simulation. The black curves shows the CMIP6 version of LMDZ for comparison with the standard 15 min time step (solid) and a shorter 1 min time step (dotted). The horizontal light gray band show the vertical ranges over which the metrics are calculated for each variable. In panel c, note that the full (resolved + subgrid) TKE from the LES is shown.</figcaption></figure>

coming from the GABLS1 LES intercomparison exercise (Beare et al., 2006), the yellow envelope in Figure 2 highlights the large range of vertical profiles obtained. This is a signature of the high sensitivity of the results to the parameters as they are varied across the range given in Table 1. In particular, very strong and unrealistic momentum decoupling manifesting as very strong wind speed gradient near the surface is allowed by the scheme in regions of the parameter space where the negative feedback of the wind shear on the mixing length (Equation 27) is overappreciated. Interestingly, such a decoupling is never simulated when using the buoyancy-only dependent length scale as shown in Figure 3b which shows the exact same experiment as in Figure 2b but using Equation 26 for the formulation of $l_s$. However, even if the yellow envelope is reasonable for the potential temperature and wind speed (Figures 3a and 3b), the use of the buoyancy-only dependent length scale can lead to unrealistically strong values of TKE in the middle of the boundary layer (Figure 3c) owing to overly high mixing length values.

Overall, the large width of the yellow envelope in Figure 2 and the possible large discrepancy with respect to the LES call for a reduction of the parameter space and a calibration of the ATKE scheme.

### 3.2.2. History Matching With Iterative Refocusing

For this purpose, we follow a history matching with iterative refocusing procedure which in practice is performed with HighTune explorer. This procedure is made of 6 steps and is fully described in Couvreux et al. (2021) and Hourdin et al. (2021). We refer the reader to the aforementioned papers for details on the method and describe here the main steps for our application.

Step 1 We first define five metrics, that is, targets for the model with respect to the LES reference, to properly capture the boundary layer structure. Those metrics are the potential temperature at the bottom (average between 30 and 60 m) and top (average between 130 and 160 m) part of the boundary layer, the zonal wind speed at the low-level jet height (average between 130 and 190 m) and the TKE at the bottom (average between 20 and 60 m) and middle (average between 60 and 100 m) part of the boundary layer. All metrics are calculated on hourly-mean profiles between the 8th and 9th hr of the simulation, when the stable boundary layer is well developed.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Same as Figure 2 but for simulations using the buoyancy length-scale formulation (Equation 26) instead of the stratification and wind-shear dependent formulation (Equation 27) in stable conditions.</figcaption></figure>

Step 2 We then define the initial parameter space consisting in a eight-dimension space corresponding to the eight parameters in bold in Figure 1 and their associated range of possible values.

Step 3 This parameter space is then sampled 80 times and experimented on GABLS1 simulation as in Section 3.2.1.

Step 4 Based on those 80 simulations, an emulator is built for each metric based on a Gaussian Process providing values for the expectation and variance at any location in the parameter space.

Step 5 We then compare the simulated metrics with respect to those from the LES reference through the calculation of an implausibility *I* for each metrics at each point *λ* of the parameter space:

$$I(\lambda) = \frac{|r - E[e_m(\lambda)]|}{\sqrt{\sigma_r^2 + \sigma_d^2 + \mathrm{Var}(e_m(\lambda))}} \quad (36)$$

where the numerator is the absolute difference between the reference metrics *r* and the corresponding expectation from the emulator *E*[*e*<sub>*m*</sub>(*λ*)]; and the denominator is the standard deviation of this difference, which includes the reference uncertainty (i.e., the spread between LES *σ*<sub>*r*</sub><sup>2</sup>), the uncertainty associated to the emulator (*Var*(*e*<sub>*m*</sub>(*λ*)), and model structural uncertainty (*σ*<sub>*d*</sub><sup>2</sup>, see Couvreux et al., 2021 for details). As the latter is not a priori known, one has to prescribe an arbitrary “tolerance to error” (see thorough discussion on the rationale behind this tolerance in Hourdin et al. (2021)) that we set to 0.25 K for potential temperature, 0.25 m s<sup>-1</sup> for wind speed and 0.01 m<sup>2</sup> s<sup>-2</sup> for TKE. History matching then rules out a part of the parameter space that corresponds to unacceptable model behavior—that is, with an implausibility higher than a given cut-off value of 3—and keeps a not-ruled out yet (NROY) space.

Step 6 Iterative refocusing then consists in sampling 80 new free parameter vectors in the NROY space and reiterates over several tuning “waves” from step 4 to 6.

Note that this procedure is not an optimization method providing in the end a single set of parameters, but a method ruling-out a non-plausible part of the initial parameter space and giving the space of acceptable free parameters—given the chosen metrics and tolerances—once it has converged.

The results after 20 waves of tuning are shown with orange envelopes for the potential temperature, wind speed and TKE profiles in Figure 2. Compared to the initial and first wave (yellow envelopes), one can first notice the convergence toward LES curves. Considerable improvement is obtained with respect to the CMIP6 version of
LMDZ, especially when the latter is run with the standard 15 min time step (solid black line) while the simulation with the 1 min time step is more reasonable. In fact, the standard LMDZ TKE-l parameterization is not by essence deficient, it seems that the numerical treatment of the TKE equation in LMDZ for this scheme does not allow convergence at time steps of about 15 min. A shallower and more realistic—compared to LES—boundary-layer height is simulated with ATKE, as well as a more peaked low-level jet and lower and much closer-to-LES TKE values. Nonetheless, the potential temperature (resp. wind speed) in the first tens of meters above the surface remains slightly overestimated (resp. underestimated). Such biases can be reduced by adding metrics targeting the lowermost part of the profiles and increasing the vertical resolution close to the surface (not shown).

We now examine the 10 “best” simulations obtained during the tuning exercise. The adjective “best” is employed here as in Hourdin et al. (2021) in the sense that the maximum (across metrics) value of the ratio of the distance to LES divided by the tolerance to error is the smallest at the end of the tuning exercise. Note that this choice of 10 simulations and the denomination “best” goes beyond the history matching philosophy as there is a priori no reason to prefer specific configurations than others in the final NROY spaces given the chosen metrics and tolerances. A choice is done here to illustrate the behavior of the ATKE scheme for single sets of parameters obtained at the end of the tuning process in 1D and 3D simulations.

Figures 4a and 4c show that they reproduce fairly well the profiles of heat and momentum turbulent fluxes, that is, two quantities that were not directly targeted during the tuning. $K_{m,h}$ values are also much lower than those in the CMIP6 physics simulation (Figures 4b and 4d) which concurs with conclusions regarding the profiles of TKE in Figure 2c. In addition, Figure 5 reveals the good numerical stability and convergence properties of the ATKE scheme in these simulations, as well as the considerable improvement regarding numerical stability and convergence with respect to the CMIP6 version of the LMDZ physics. One can in fact point out the serious problems of numerical convergence and stability of the TKE that was present in the CMIP6 version of the LMDZ physics in weakly stable boundary layer conditions. In fact it appears that a substantial part of the better performances of the ATKE scheme compared to the CMIP6 version of the LMDZ physics in Figures 2 and 4 come from the better numerical properties of ATKE. This analysis overall makes us confident with the robustness and efficiency of the numerical resolution method for the TKE evolution equation in the ATKE scheme presented in Section 2.2.2.

Figure 6 shows the implausibility matrices after 20 waves of history matching. We refer the reader to Hourdin et al. (2021) for an exhaustive description of this type of figure. The matrix is divided into 2D sub-matrices, each one being a restriction to 2 parameters, the names of which are given along the main diagonal. Each axis spans the initial [min,max] ranges for the $n = 8$ parameter considered and shown in Table 1. The sub-matrices of the lower-left triangle show for each pixel the minimum implausibility obtained when varying the $n-2$ other parameters. The matrices in upper right triangles show fraction of points with implausibility lower than the cut-off value, when varying the $n-2$ other parameters. Gray color means that the reference LES cannot be matched by the simulations by varying the $n-2$ unfixed parameters while a value of 1 (yellow color) means that values of the two parameters in x and y axis can be retained whatever the values of the $n-2$ other parameters. When inspecting the NROY space after 20 waves of tuning, one can notice that its final shape has been mostly constrained by the $c_l$ and $c_e$—only high values of those two parameters are kept in the NROY space—and to a lesser extent by $l_\infty$. This does not absolutely mean that the other five parameters do not play role in the overall behavior of the scheme but this shows that the representation of the GABLS1 weakly stable boundary layer with ATKE mostly depends upon the value of $c_l$, $c_e$, and $l_\infty$. This point is further shown by the strong similarity between Figure 7—which has been produced with a tuning on $c_l$, $c_e$ and $l_\infty$ only—and Figure 2. Such a result is not that surprising since the analysis of the terms of the TKE (not shown) reveals that the TKE during GABLS1 results from a quasi-equilibrium between wind shear production and dissipation. Therefore, the simulation of the turbulent diffusion in this weakly stable boundary layer necessarily mostly depends on the parameters that control the magnitude of the mixing length and the intensity of the shear production and TKE dissipation namely $c_l$, $c_e$, and $l_\infty$.

It is worth mentioning that we have carried out another tuning exercise (not shown here) which also considered metrics on typical cloudy and convective boundary layers 1D cases as well as a more extensive sets of parameters including parameters from the cloud and mass-flux schemes of LMDZ, in line with Section 5.1 of Hourdin et al. (2021). Interestingly, this tuning also evidenced the important role of $c_l$ and $c_e$ and the History Matching procedures kept only high values for those two parameters (similarly to Figure 6). This result suggests that $c_l$ and
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Vertical profiles of momentum flux (panel a), heat flux (panel c), eddy diffusivity coefficient for momentum (panel b) and heat (panel d) after 9 hr of GEWEX Atmospheric Boundary Layer Study 1 simulation. Gray curves show the LMDZ simulations run with the 10 best parameter vectors after the tuning exercise. Blue curves in panels (a) and (c) show the five reference Large Eddy Simulations. The red curve shows the “best” LMDZ simulation obtained during the tuning exercise (see main text for details). The solid (resp. dotted) black line shows the CMIP6 version of LMDZ with the standard 15 min time step (resp. with a shorter 1 min step) for comparison.</figcaption></figure>

$c_e$ are mostly constrained by the metrics on the stable boundary layers than those on convective boundary layers for which or which the mixing is mostly controlled by the LMDZ mass-flux scheme.

The overall weak dependence upon $c_e$ shown in Figure 6 may have somewhat been expected given the relatively weak contribution of the transport term $\mathcal{T}$ is the overall TKE budget (not shown). Regarding $S_{\text{min}}$, $Ri_c$, and $\alpha_{\text{Pr}}$, one may expect a more important role of those parameters in very stable boundary layers that is, with a stratification
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Time evolution of the Turbulent Kinetic Energy at 40 m a.g.l. in LMDZ single column model GEWEX Atmospheric Boundary Layer Study 1 simulations. Solid gray curves show the simulations run with the 10 best parameter vectors after the tuning exercise and a 15 min time step. The solid and dotted red curves shows simulations run with the best parameter vector and a time step of 15 and 1 min respectively. The solid and dotted black curves shows simulations run with CMIP6 version of LMDZ and a time step of 15 and 1 min respectively.</figcaption></figure>

more pronounced compared to that in GABLS1. Their values might thus be more constrained if we were to tune the ATKE scheme over a more stable boundary layer case such as GABLS4 (Couvreux et al., 2020) instead of or in addition to GABLS1. However LES do not converge that well on GABLS4 which makes the tuning exercise more delicate. Moreover, the role of radiation in determining the structure of the boundary-layer becomes increasingly important as stability increases (Edwards, 2009) and in addition to turbulent diffusion, the coupling between turbulence and radiation becomes an essential feature to capture with models. We therefore leave this aspect for further research.

### 3.3. Challenging the Antarctic and Martian Stable Boundary Layers

We now conduct two short and arbitrary applications of the ATKE parameterization in simulations with the LMDZ GCM and Mars PCM.

#### 3.3.1. Stable Boundary Layer Regimes at Dome C, Antarctic Plateau

First, we verify that the proposed scheme is able to reproduce the dichotomous behavior of the stable boundary layer at Dome C on the Antarctic Plateau that is, a very stable regime with strong temperature surface-based inversions and collapsed turbulence versus a weakly stable state with weak inversions. The sharp transition between those two regimes occurs in a narrow range of wind speed (Baas et al., 2019; Vignon, van de Wiel, et al., 2017) and can trigger the formation of ice fog clouds (Vignon et al., 2022). Such a test was proposed in Vignon et al. (2018) to verify the ability of the CMIP6 version of LMDZ to reproduce the overall dynamics of the stable boundary layers and it is performed here as capturing the Dome C boundary layer was identified as a target during the development of LMDZ for CMIP6 (Cheruy et al., 2020). This is an aspect that we want to conserve throughout the development of the LMDZ physics and particularly when introducing a new turbulent diffusion scheme. It is also worth noting that such a test was also used for the recent development of the CanAM model (He et al., 2019) as well as for verifying the robustness of LES of the stable boundary layer (van der Linden et al., 2019). We follow here the exact same LMDZ simulation configuration as in Vignon et al. (2018) that is, 1 year (2015) simulations are conducted with the zooming capability of the LMDZ to refine a 64 × 64 global grid to reach a 50 × 50 km on the Dome C. One slight difference though with respect to Vignon et al. (2018) is that we use the 95-level vertical grid used in the previous section instead of the 79-level grid in the reference paper. It is worth mentioning that sensitivity experiments on GABLS1 reveal that once properly and specifically tuned for the vertical grid considered, the performances of the ATKE scheme in simulating the stable boundary layer is quite similar when using either the 79-level or the 95-level vertical grid. Nudging in wind, temperature and humidity toward ERA5 reanalyzes (Hersbach et al., 2020) is applied outside the zoom area to evaluate the sub-components of the physics of the model apart from likely deficiencies in representing the large scale meteorological fields. The reader is referred to Vignon et al. (2018) for details on the simulation configuration and the surface snow treatment in LMDZ and to Genthon et al. (2013) and Vignon, van de Wiel, et al. (2017) for information about the vertical structure of the boundary layer at Dome C. The simulation has been run with the CMIP6 version of the LMDZ physics as well as by an adapted versions using the ATKE diffusion scheme and the 10 "best" sets of parameters found from the single column model tuning.

A simple diagnostics to assess the representation of the two stable boundary layer regimes is to investigate the dependence of the surface-based temperature inversion upon the wind speed in clear sky conditions. Data align along a well-defined "inverted-S" shape curve (van de Wiel et al., 2017; Vignon, van de Wiel, et al., 2017), the two horizontal branches corresponding to the two regimes and the vertical one to the non-linear transition between them as the wind speed increases or decreases (Figure 8a). In the very stable regime, the surface-layer Richardson number can reach values up to 20. As shown in Figure 8b, the CMIP6 version of LMDZ reasonably captures the strong surface-atmosphere decoupling in very stable conditions and the two-regime behavior. LMDZ with the ATKE scheme run with the "best" set of parameters (Figure 8c) retained in Section 3.2 reproduces even more realistically reproduce the two-regime behavior—that is, the reversed "S" shape pattern—and the decoupling in
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Implausibility matrix after 20 waves of history matching exploration. The upper-right triangle is made of sub-matrices that show the fraction of points with implausibility lower than the chosen cutoff while the sub-matrices of the lower-left triangle show the minimum value of the implausibility when all the parameters are varied except those used as x- and y-axis, the name of which are given on the diagonal of the main matrix. The number at the bottom of the graph shows the not-ruled out yet space value (fraction of the initial parameter space) after 20 waves.</figcaption></figure>

very stable conditions despite an overestimation of the strong near-surface temperature inversions. The latter can be attributed to an overly weak downward longwave radiative flux from the very dry and cold Dome C atmosphere in clear-sky conditions (Vignon et al., 2018). Not only the two-regime transition and the “inverted S” shape is reproduced by the LMDZ with ATKE but the mean vertical profiles of temperature inversion and wind speed in both regimes are also reasonably well simulated (Figure 9). Persistent biases are the slightly underestimated wind speed in the weakly stable regime (Figure 9c) and higher temperature values—with respect to the ground surface—compared to observations above 20 m in the very stable regime (Figure 9b), which concurs with the overestimation of the strong near-surface temperature inversions.

An important point here is that the two-regime behavior is captured with all the 10 “best” sets of parameters after 20 waves of tuning on GABLS1 (Figures 8c–8l) and despite the fact that such a GABLS1-based tuning has not substantially constrained parameters that may be a priori important in very stable conditions such as $S_{min}$, $Ri_c$, and $\alpha_{Pr}$. In fact, the transition between the weakly and very stable regimes of the stable boundary-layer primarily relies on the ability of a TKE-l scheme to allow for a turbulence collapse in very stable conditions (Vignon et al., 2018).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Figure 2 but after a tuning on $c_e$, $c_l$ and $l_\infty$ only. The other parameters have been arbitrarily set to the following values: $Ri_c = 0.2$, $S_{min} = 0.05$, $Pr_n = 0.8$, $\alpha_{Pr} = 4.5$ and $c_e = 2.0$. Note that we have stopped the tuning exercise at the 9th wave here since convergence has been attained.</figcaption></figure>

This is the case with the ATKE scheme—whatever the $S_{min}$, $Ri_c$, and $\alpha_{Pr}$ value chosen in their corresponding ranges of acceptable values—as no artificial threshold or lower-bound has been prescribed to maintain a certain amount of TKE in very stable conditions.

### 3.3.2. Nocturnal Stable Boundary Layer Collapse on Mars

Mars has a thinner and much less dense atmosphere compared to Earth and its planetary boundary layer exhibits stronger diurnal variations (Petrosyan et al., 2011; Spiga et al., 2010) with a abrupt collapse at the day-night transition. During night-time, the Martian boundary layer exhibits numerous similarities with that of the polar regions on Earth such as strong surface-based temperature inversions associated with very weak turbulence (Banfield et al., 2020), the latter being able to re-activate through wind shear production associated with low-level jets (Chatain et al., 2021).

This extreme environment enables us to challenge the versatility of ATKE parameterization and compare its performance with the default TKE-1 scheme used in the current Mars PCM (Colaïtis et al., 2013).

As a first test, we compare the two parameterizations using the single-column version of the Mars PCM to assess the overall behavior of the diurnal cycle of the boundary layer and the numerical stability of the model. The single-column version of the Mars PCM uses the same physics as the 3D model (Lange et al., 2023) in which the ATKE scheme has been implemented as an option and a vertical grid with six levels in the first km above the ground. No lateral advection of heat and momentum is prescribed, the initial temperature profile is set to 180 K and the zonal wind speed is nudged toward a constant value of 7 m s⁻¹ which corresponds to values measured at the Mars Equator by the InSight lander (Banfield et al., 2020). Simulations are performed at the Equator, with no dust aerosols, and ran for several Martian days until the diurnal cycle reaches an equilibrium after 10 days. The nocturnal boundary layer simulated is weakly to moderately stable, with a near-surface gradient Richardson reaching a value of ≈ 0.1. Figure 10 shows the evolution of the TKE (color shading) and wind speed (contours) in the first km above the ground surface during a typical diurnal cycle. As explained in Section 3.1, the nocturnal TKE field simulated by the default TKE-1 scheme of the Mars PCM is affected by strong numerical oscillations (Figure 10a) which are mitigated when adding a minimum mixing coefficient $K_{min}$ (Figure 10b). When running the Mars PCM using the ATKE scheme with the “best” set of parameters retained from the tuning on GABLS1 in Section 3.2.2 (Figure 10c) and with no prescription of $K_{min}$, the structure of the nocturnal boundary layer is well captured and no numerical oscillations affect the TKE and wind fields. Unlike in Figure 10b, the TKE exhibits a
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Dome C temperature inversion between 10 m and the ground surface plotted as a function of the 10-m wind speed in clear-sky conditions (downward longwave radiative flux <100 W m⁻²) from April to September 2015. Panel (a) shows results from in situ observations. Panel (b) (resp. c) show the LMDZ simulation in the CMIP6 physics configuration (resp. with the ATKE scheme using the “best” set of parameters retained in Section 3.2). Panels (d) to (l) show results from nine simulations with the ATKE scheme using nine following “best” sets of parameters after the tuning phase on GEWEX Atmospheric Boundary Layer Study 1. Dome C measurement data are from Genthon et al. (2021).</figcaption></figure>
continuous decrease with increasing height in the nocturnal boundary layer, which better concurs with the typical TKE structure in weakly stable boundary layers (e.g., Acevedo et al., 2015).

We then assess the performance of the ATKE model by performing simulations with the 3D Mars PCM and comparing the results to in situ wind observations collected by the InSight lander deployed at a latitude 4.5°N and a longitude of 135°E. InSight continuously monitored the wind at a height of 1.2 m for almost one martian year with an unprecedented time resolution (Banfield et al., 2020). Two striking phenomena have been detected. First, a dramatic reduction of the wind speed, following the collapse of the boundary layer is observed around 17–18 local time during the clear season (Figure 11a) that is, the first half of the Martian year when a relatively small amount of dust is present in the Martian sky (Kahre et al., 2017). The abruptness of this change is related to both the very low thermal inertia of the Martian ground surface and the thinness of the Martian atmosphere. Second, during the dusty season that is, the second half of the Martian year, substantial night-time turbulence is observed (Chatain et al., 2021) and the decrease in near-surface wind speed is less pronounced (Figure 11d). Those two phenomena have been shown to be poorly reproduced by the Mars PCM, in particular, the collapse of winds at sunset (Forget et al., 2021).

Here, as a proof of concept, we run the 3D Mars PCM using either the default TKE-1 scheme and the ATKE scheme with the “best” set of parameters from the GABLS1 tuning that is, with no specific tuning for Martian conditions. We run the 3D global model over one complete martian year and with a horizontal resolution of 3.75° in latitude and 5.625° in longitude. Initial conditions are derived from 10-year simulations which provide equilibrium states of water and CO₂ cycles (Pottier et al., 2017). The seasonal and geographic variations of dust opacity in the sky are prescribed using dust observations by Montabone et al. (2015). Results are presented in Figure 11. Concurring with Forget et al. (2021), the model in its standard configuration fails to reproduce the sharp transition from high to low wind speeds at sunset during clear conditions (Figure 11b). This aspect is significantly improved when using the ATKE scheme (Figure 11c). However, the wind speed in the second part of the night remains underestimated in both configurations which questions the representation of the surface-atmosphere coupling in this period (Chatain et al., 2021). In particular, the ATKE scheme produces a close-to-zero mean wind speed a few hours after sunset. Further investigation reveals that the overall wind speed in Mars' troposphere weakens—meaning that the large-scale pressure force is weak—and the near-surface Richardson number becomes sufficiently high such that almost no TKE is maintained near the surface leading to an overly weak top-down transfer of momentum, a common mechanism in very stable boundary layers (Acevedo et al., 2015). In fact, Chatain et al. (2021) show that the turbulence during this period of the night is mostly non-local and results from meso-scale structures which cannot be reproduced by the current version of the ATKE scheme.

In the dusty season, the current model overestimates the surface wind speed owing to an excess of turbulent mixing (Figure 11e), while the ATKE parameterization leads to more realistic wind speeds (Figure 11f).

Overall, this preliminary experiment demonstrates: (a) the applicability of the ATKE parameterization on Mars and the promising results that can be obtained with a set of parameters not specifically tuned for Mars conditions and; (b) the improvement of the model both numerically and physically in stable conditions. Nonetheless, Mars simulations with the ATKE scheme would further benefit from a more adapted tuning using references such as Mars LES (Spiga et al., 2010) or InSight observations (Banfield et al., 2020). It is also worth noting that the Mars atmosphere, particularly at the poles that is, far from the InSight landing site, exhibits particularities that cannot be properly captured with the current version of the ATKE scheme. A key aspect is that air buoyancy can be created by compositional vertical gradients of both water vapor and carbon dioxide, that is, the prevailing gas of Mars' atmosphere. In particular, during the winter polar night, CO₂ condenses upon the ice cap surface (e.g., Weiss & Ingersoll, 2000) changing dramatically the near-surface atmospheric composition. Such an effect cannot be taken into account given with Brünt-Vaisala pulsation and Richardson number expressions based on a virtual potential temperature. This aspect deserves attention for further improvement of the ATKE scheme.

## 4. Summary and Conclusions

This study presents the development of a simple TKE-1 parameterization of turbulent eddy coefficients for the simulation of the neutral and stable boundary layer in large-scale atmospheric models. The parameterization has been carefully designed such that all adjustable parameters have been clearly identified and their ranges of possible values defined to help the calibration and assess the parametric sensitivity. Instead of using fixed and
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Dome C mean vertical profiles of surface-based temperature inversion (left) and wind speed (bottom row) in the weakly stable (left) and very stable regimes of the stable boundary layer in clear sky conditions. Weakly (resp. very) stable situations which correspond to the lower (resp. higher) horizontal branches of the “inverted S” in Figure 8 are defined when the difference in temperature between 10 m and the surface is lower than 5 K (resp. higher than 20 K). Orange curves show the in-situ observations, the black curves the LMDZ simulation in the CMIP6 physics configuration, the red curve the LMDZ with the ATKE scheme using the “best” set of parameters retained in Section 3.2 and the gray lines the results from nine simulations with the ATKE scheme using nine following “best” sets of parameters after the tuning phase on GEWEX Atmospheric Boundary Layer Study 1.</figcaption></figure>

empirical expressions of stability functions and turbulent Prandlt number, we have derived fully tunable and heuristic formulae to improve the versatility of the scheme and its potential applicability for planetary atmospheres composed of an ideal and perfect gas. A wind-shear and buoyancy dependent formulation for the mixing length in stratified conditions is considered. A two-step numerical treatment of the TKE equation is further proposed and shows good convergence and stability properties at typical time steps used in large scale atmospheric models. Despite its first order importance to simulate the stable boundary layers in atmospheric models, this aspect has not been really addressed in the well-known stable boundary layer intercomparison exercises (Bazile et al., 2014; Bosveld et al., 2014; Cuxart et al., 2006; Svensson et al., 2011) and probably deserves more attention. The parametric sensitivity of the ATKE scheme has been assessed with the HighTune explorer tools using 1D simulations of the GABLS1 weakly stable boundary layer case with the single-column version of LMDZ. Using a History-Matching approach, we carried out a first calibration of the scheme allowing us to reduce the initial parameter space to keep an ensemble that satisfies the representation of weakly stable boundary layers.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Evolution of the Turbulent Kinetic Energy through the Martian day in (a) the baseline physics configuration; (b) the same configuration with the inclusion of a minimum mixing coefficient $K_{min}$; (c) the simulation using the ATKE scheme for turbulent diffusion. Black contours indicate the wind speed in m s⁻¹.</figcaption></figure>

Substantial improvement with respect to the CMIP6 version of LMDZ has been achieved in terms of vertical profiles of temperature, wind, TKE and turbulent fluxes of momentum and heat. Such improvements are first explained by the better numerical treatment of the TKE equation and to a second extent, by the calibration of parameters—especially those controlling the value of the mixing length, TKE shear production and dissipation

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison between InSight wind speed measurements (gray dots and black curves in panels a and d) and Mars PCM simulations using the default TKE-I scheme (b, e) and the ATKE scheme (c, f). For model fields, the mean wind speed over the period considered is presented in solid lines, and the diurnal variability is shown with the envelope of dashed lines ($q_1$ and $q_3$ referring to the first and third quartiles).</figcaption></figure>
terms—which can be extensively performed with ATKE thanks to the chosen formulations of the mixing length, stability functions and turbulent Prandtl number.

However the tuning experiment restricted to the weakly stable GABLS1 case has not enabled us to clearly evidence a potential added value of a wind-shear and buoyancy dependent formulation for the mixing length in stratified conditions compared to a buoyancy only-dependent one. The ability of the ATKE scheme to simulate the stable boundary layer as well as its applicability to planetary atmospheres have then been assessed through simulations of the Antarctic and Martian boundary layer with the LMDZ and Mars Planetary Climate model respectively. In particular the two-regime behavior of the stable boundary layer at Dome C, a challenge for turbulent diffusion schemes in GCMs, is reasonably well captured with the ATKE scheme. In addition, promising results have been obtained for the representation of the nocturnal Martian boundary layer with improvements regarding the numerical stability compared to the original model. Such results pave the way for a Mars-specific tuning of the ATKE scheme in the future.

A prospect of our work is to verify the physical and numerical robustness of the ATKE parameterization in atmospheric flows with extremely strong wind shear such as katabatic winds developing over ice caps. Such an application could also make it possible to assess a potential added value of a wind shear-dependent formulation of the mixing length. Moreover, in view of a fully reliable application in a climate model such as LMDZ, the key parameters of the ATKE scheme—especially $c_l$ and $c_e$—should be included in a more thorough tuning exercise including parameters from other parameterizations and considering additional metrics on convective boundary layer simulations (Hourdin et al., 2021). It is worth mentioning that additional and exploratory 1D test simulations of typical convective boundary layer cases with the ATKE parameterization—with the “best” set of parameters after the tuning exercise on GABLS1—coupled with the LMDZ mass-flux scheme show promising results and similar performances with respect to the CMIP6 physics of LMDZ (not shown here).

Last but not least, we would like to emphasize that this work was initiated and fostered during collaborative work sessions dedicated to the transfer of knowledge and critical questioning on the physics and assumptions behind the parameterizations used in planetary GCMs. Those sessions spontaneously emerged following students' questions and gathered atmospheric and planetary scientists experts and non experts of turbulent mixing and parameterization development. The motivations behind the ATKE scheme development went beyond the need to advance the turbulent diffusion scheme in our models but were also—and maybe first—a reason and a need to teach and learn the parameterization development in a “learning-by-doing” way.

## Appendix A: A Gravity-Invariant Formulation of Our TKE-l Turbulent Diffusion Scheme

For the sake of universality of a turbulent diffusion parameterization and in particular for potential application on different planets, one may want to develop a framework as independent as possible upon planet's characteristics, in particular upon planet's gravity. In the main paper, gravity appears in the expression of the Brünt Väisälä frequency thus in the expression of the gradient Richardson number and in the buoyancy term of the Turbulent Kinetic Energy (TKE) evolution equation Equation 7. In this appendix, we briefly introduce a framework using geopotential as vertical coordinate and in which gravity is no longer involved. Such a framework is proposed here as a prospect for a further new implementation of the parameterization.

Let's introduce the geopotential $\phi$ defined such that $d\phi = gdz$ as well as a “re-scaled” time $\tau$ defined by $d\tau = gdt$. The diffusion equation of a quantity $c$ (Equation 5) can be written in the form:

$$ \frac{\partial c}{\partial \tau} = \frac{1}{\rho} \frac{\partial}{\partial \phi} \left( \rho K_c^{\phi} \frac{\partial c}{\partial \phi} \right) \qquad (A1) $$

where $K_c^{\phi} = gK_c$. In such a framework, assuming down-gradient expression of turbulent fluxes and the same closures for the TKE dissipation and transport terms as in the main manuscript, the TKE evolution Equation A1 reads:

$$ \frac{\partial e}{\partial \tau} = K_m^{\phi} \left[ (S^{\phi})^2 - \text{Pr}(Ri)(N^{\phi})^2 \right] + \frac{1}{\rho} \frac{\partial}{\partial \phi} \left( \rho c_e K_m^{\phi} \frac{\partial e}{\partial \phi} \right) - \frac{e^{3/2}}{c_e l^{\phi}} \qquad (A2) $$
with $l^\phi = gl$, $(S^\phi)^2 = (\partial_\phi u)^2 + (\partial_\phi v)^2$, and $(N^\phi)^2 = \frac{1}{\partial_\phi} \frac{\partial \theta_\phi}{\partial \phi}$.

One can then express $K_m^\phi = l^\phi(\phi, e, Ri)S_m(Ri)\sqrt{e}$. Noting the gravity independent form of the gradient Richardson number $Ri = (N^\phi)^2/(S^\phi)^2$, the expressions for $S_m(Ri)$ and $Pr(Ri)$ can be taken identically from Equations 20 and 23 as they are gravity-independent. For the mixing length $l^\phi$ expression, one can use a similar approach as in Section 2.4 replacing the neutral-limit formulation with

$$l_n^\phi = \frac{\kappa \phi l_\infty^\phi}{\kappa \phi + l_\infty^\phi} \qquad (A3)$$

$l_\infty^\phi$ being a tuning parameter. In such a way, Equations A1 and A2 combined with the proposed expressions for $K_m$, Pr and $l^\phi$ establish a complete gravity-invariant formulation of the turbulent diffusion parameterization.

## Data Availability Statement

The latest version of the LMDZ source code can be downloaded freely from the LMDZ web site. The version used for the specific simulation runs for this paper is the “svn” release 4781 from 21 December 2023, which can be downloaded and installed on a Linux computer by running the install_lmdz.sh script available here: https://web.lmd.jussieu.fr/~lmdz/pub/pub/install_lmdz.sh. The Mars PCM used in this work can be downloaded with documentation from the SVN repository at https://svn.lmd.jussieu.fr/Planeto/trunk/LMDZ.MARS/. Forcings for the GABLS1 single-column cases are provided under the DEPHY-SCM standard at the following link: https://github.com/GdR-DEPHY/DEPHY-SCM/. GABLS1 LES used in the intercomparison exercise of Beare et al. (2006) are distributed here: https://gabls.metoffice.gov.uk/lem_data.html.

Dome C temperature and wind speed data are freely distributed in Genthon et al. (2021). InSight wind data can be retrieved from the Planetary Data System (Jose Rodriguez-Manfredi, 2019).

## Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant 951596) through the AWACA project. The part of the work related to the Martian atmosphere is supported by the ANR project MAGIS. We acknowledge support from the DEPHY research group, funded by CNRS/INSU and Météo-France, as well as from the PEPR TRACCS project (no. ANR-22-EXTR-0008) funded from the Agence Nationale de la Recherche—France 2030. Éric Bazile is gratefully acknowledged for fruitful discussions. This study was conducted using the ESPRI (Ensemble de Services Pour la Recherche l'IPSL) computing and data center (https://mesocentre.ipsl.fr) which is supported by CNRS, Sorbonne Université, Ecole Polytechnique, and CNES and through national and international grants. Simulations were performed using HPC resources from the IDRIS (Institut du Développement et des Ressources en Informatique Scientifique, CNRS, France), project RLMD AD010107632R1. Simulations with the Mars PCM were performed using HPC resources of Centre Informatique National de l'Enseignement Supérieur (CINES) under the allocation no A0100110391 made by Grand Equipment National de Calcul Intensif (GENCI). We gratefully thank three anonymous reviewers for their careful evaluation of the manuscript and thoughtful comments.

## References

* Acevedo, O. C., Mahrt, L., Puhales, F. S., Costa, F. D., Medeiros, L. E., & Degrazia, G. A. (2015). Contrasting structures between the decoupled and coupled states of the stable boundary layer. *Quarterly Journal of the Royal Meteorological Society*, 142(695), 693–702. https://doi.org/10.1002/qj.2693

* André, J. C., Moor, G. D., Lacarrère, P., & du Vachat, R. (1978). Modeling the 24-hour evolution of the mean and turbulent structures of the planetary boundary layer. *Journal of the Atmospheric Sciences*, 35(10), 1861–1883. https://doi.org/10.1175/1520-0469(1978)035(1861:MTHEOT)2.0.CO;2

* Audouin, O., Roehrig, R., Couvreux, F., & Williamson, D. (2021). Modeling the gabls4 strongly-stable boundary layer with a gcm turbulence parameterization: Parametric sensitivity or intrinsic limits? *Journal of Advances in Modeling Earth Systems*, 13(3), e2020MS002269. https://doi.org/10.1029/2020MS002269

* Baas, P., Van De Wiel, B., Van der Linden, S., & Bosveld, F. (2018). From near-neutral to strongly stratified: Adequately modelling the clear-sky nocturnal boundary layer at cabauw. *Boundary-Layer Meteorology*, 166(2), 217–238. https://doi.org/10.1007/s10546-017-0304-8

* Baas, P., van de Wiel, B. J. H., van Meijgaard, E., Vignon, E., Genthon, C., van der Linden, S. J. A., & de Roode, S. R. (2019). Transitions in the wintertime near-surface temperature inversion at dome c, Antarctica. *Quarterly Journal of the Royal Meteorological Society*, 145(720), 930–946. https://doi.org/10.1002/qj.3450

* Banfield, D., Spiga, A., Newman, C., Forget, F., Lemmon, M., Lorenz, R., et al. (2020). The atmosphere of Mars as observed by insight. *Nature Geoscience*, 13(3), 190–198. https://doi.org/10.1038/s41561-020-0534-0

* Bazile, E., Couvreux, F., Le Moigne, P., Genthon, C., Holtslag, A. A. M., & Svensson, G. (2014). GABLS4: An intercomparison case to study the stable boundary layer over the Antarctic Plateau. *GEWEX News*, 24, 4.

* Bazile, E., Marquet, P., Bouteloup, Y., & Bouyssel, F. (2011). The turbulent kinetic energy (TKE) scheme in the nwp models at météo-France. In *Workshop on diurnal cycles and the stable boundary layer, ECMWF* (pp. 127–136).

* Beare, R. J., Macvean, M. K., Holtslag, A. A. M., Cuxart, J., Esau, I., Golaz, J.-C., et al. (2006). An intercomparison of large-eddy simulations of the stable boundary layer. *Boundary-Layer Meteorology*, 118(2), 247–272. https://doi.org/10.1007/s10546-004-2820-6

* Betts, A. K. (1973). Non-precipitating cumulus convection and its parameterization. *Quarterly Journal of the Royal Meteorological Society*, 99(419), 178–196. https://doi.org/10.1002/qj.49709941915

* Blackadar, A. K. (1962). The vertical distribution of wind and turbulent exchange in neutral atmosphere. *Journal of Geophysical Research*, 67(8), 3095–3102. https://doi.org/10.1029/jz067i008p03095

* Bosveld, F. C., Baas, P., van Meijgaard, E., de Bruijn, E. I. F., Steeneveld, G.-J., & Holtslag, A. A. M. (2014). The third GABLS intercomparison case for evaluation studies of boundary-layer models. Part A: Case selection and set-up. *Boundary-Layer Meteorology*, 152(2), 133–156. https://doi.org/10.1007/s10546-014-9917-3

* Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., et al. (2020). Presentation and evaluation of the IPSL-CM6A-LR climate model. *Journal of Advances in Modeling Earth Systems*, 12(7), e2019MS002010. https://doi.org/10.1029/2019MS002010

* Bougeault, P., & Lacarrère, P. (1989). Parametrization of orography-induced turbulence in a mesobeta-scale model. *Monthly-Weather Reviews*, 117, 1872–1891.
Chatain, A., Spiga, A., Banfield, D., Forget, F., & Murdoch, N. (2021). Seasonal variability of the daytime and nighttime atmospheric turbulence experienced by insight on Mars. *Geophysical Research Letters, 48*(22), e2021GL095453. https://doi.org/10.1029/2021GL095453

Chen, W., Lovejoy, S., & Muller, J.-P. (2016). Mars' atmosphere: The sister planet, our statistical twin. *Journal of Geophysical Research: Atmospheres, 121*(20), 11968–11988. https://doi.org/10.1002/2016JD025211

Cheruy, F., Ducharne, A., Hourdin, F., Musat, I., Vignon, E., Gastineau, G., et al. (2020). Improved near-surface continental climate in IPSL-CM6A-LR by combined evolutions of atmospheric and land surface physics. *Journal of Advances in Modeling Earth Systems, 12*(10), e2019MS002005. https://doi.org/10.1029/2019MS002005

Colaïtis, A., Spiga, A., Hourdin, F., Rio, C., Forget, F., & Millour, E. (2013). A thermal plume model for the Martian convective boundary layer. *Journal of Geophysical Research: Planets, 118*(7), 1468–1487. https://doi.org/10.1002/jgre.20104

Couvreux, F., Bazile, E., Rodier, Q., Maronga, B., Matheou, G., Chinita, M. J., et al. (2020). Intercomparison of large-eddy simulations of the Antarctic boundary layer for very stable stratification. *Boundary-Layer Meteorology, 176*(3), 369–400. https://doi.org/10.1007/s10546-020-00539-4

Couvreux, F., Hourdin, F., Williamson, D., Roehrig, R., Volodina, V., Villefranque, N., et al. (2021). Process-based climate model development harnessing machine learning: I. A calibration tool for parameterization improvement. *Journal of Advances in Modeling Earth Systems, 13*(3), e2020MS002217. https://doi.org/10.1029/2020MS002217

Cuxart, J., Holtslag, A. A. M., Beare, R. J., Bazile, E., Beljaars, A., Cheng, A., et al. (2006). Single-column model intercomparison for a stably stratified atmospheric boundary layer. *Boundary-Layer Meteorology, 118*(2), 273–303. https://doi.org/10.1007/s10546-005-3780-1

Deardoff, J. W. (1980). Stratocumulus-capped mixed layers derived from a three dimensional model. *Boundary-Layer Meteorology, 18*(4), 495–527. https://doi.org/10.1007/bf00119502

Delage, Y. (1997). Parametrising sub-grid scale vertical transport in atmospheric models under statically stable conditions. *Boundary-Layer Meteorology, 82*(1), 23–48. https://doi.org/10.1023/a:1000132524077

Deleersnijder, E. (1992). *Modélisation hydrodynamique tridimensionnelle de la circulation générale estivale de la région du détroit de bering (in french)* (Unpublished doctoral dissertation). Université Catholique de Louvain.

Derbyshire, S. H. (1990). Nieuwstadt's stable boundary layer revisited. *Quarterly Journal of the Royal Meteorological Society, 126*(491), 127–158. https://doi.org/10.1256/smsqi.49105

Dufresne, J.-L., & Ghattas, J. (2009). Description du schéma de la couche limite turbulente et l'interface avec la surface planetaire dans lmdz. In *ORCHIDEE documentation, in French*. Retrieved from http://forge.ipsj.lussieu.fr/orchide/attachment/wiki/Documentation/CouplingLMDZ/Dufresne,%20Ghattas%20-%202009_Coupling-ORC-LMDZ.pdf

Edwards, J. M. (2009). Radiative processes in the stable boundary layer: Part i. radiatives aspects. *Boundary-Layer Meteorology, 131*(1), 105–126. https://doi.org/10.1007/s10546-009-9364-8

Emanuel, K. A. (1991). A scheme for representing cumulus convection in large-scale models. *Journal of the Atmospheric Sciences, 48*(21), 2313–2329. https://doi.org/10.1175/1520-0469(1991)048(2313:ASFRCC)2.0.CO;2

England, D. E., & McNider, R. T. (1995). Stability functions based upon shear functions. *Boundary-Layer Meteorology, 74*(1–2), 113–130. https://doi.org/10.1007/bf00715713

Forget, F., Banfield, D., Spiga, A., Millour, E., Borella, A., Lange, L., et al. (2021). More than one Martian year of meteorology observed by the InSight Lander. In *European planetary science congress EPSC2021-273*. https://doi.org/10.5194/epsc2021-273

Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., et al. (1999). Improved general circulation models of the martian atmosphere from the surface to above 80 km. *Journal of Geophysical Research, 104*(E10), 24155–24175. https://doi.org/10.1029/1999JE001025

Genthon, C., Six, D., Gallée, H., Grigioni, P., & Pellegrini, A. (2013). Two years of atmospheric boundary layer observations on a 45-m tower at Dome C on the Antarctic Plateau. *Journal of Geophysical Research: Atmospheres, 118*(8), 3218–3232. https://doi.org/10.1029/jgrd.50128

Genthon, C., Veron, D., Vignon, E., Six, D., Dufresne, J.-L., Madeleine, J.-B., et al. (2021). 10 years of temperature and wind observation on a 45 m tower at Dome C, East Antarctic plateau. *Earth System Science Data, 13*(12), 5731–5746. https://doi.org/10.5194/essd-13-5731-2021

Girard, C., & Delage, Y. (1990). Stable schemes for nonlinear vertical diffusion in atmospheric circulation models. *Monthly Weather Review, 118*(3), 737–745. https://doi.org/10.1175/1520-0493(1990)118(0737:SSFNVVD)2.0.CO;2

Golaz, J.-C., Larson, V. E., & Cotton, W. R. (2002). A PDF-based model for boundary layer clouds. Part I: Method and model description. *Journal of the Atmospheric Sciences, 59*(24), 3540–3551. https://doi.org/10.1175/1520-0469(2002)059(3540:APBMFB)2.0.CO;2

Grisogono, B. (2010). Generalizing 'z-less' mixing length for stable boundary layers. *Quarterly Journal of the Royal Meteorological Society, 136*(646), 213–221. https://doi.org/10.1002/qj.529

Grisogono, B., & Belušić, D. (2008). Improving mixing length-scale for stable boundary layers. *Quarterly Journal of the Royal Meteorological Society, 134*(637), 2185–2192. https://doi.org/10.1002/qj.347

He, Y., McFarlane, N. A., & Monahan, A. H. (2019). A new tke-based parameterization of atmospheric turbulence in the Canadian global and regional climate models. *Journal of Advances in Modeling Earth Systems, 11*(5), 1153–1188. https://doi.org/10.1029/2018MS001532

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society, 146*(730), 1999–2049. https://doi.org/10.1002/qj.3803

Holtslag, A. A. M., & Boville, B. A. (1993). Local versus non-local boundary layer diffusion in a global climate model. *Journal of Climate, 6*(10), 1825–1842. https://doi.org/10.1175/1520-0442(1993)006<1825:lvnbld>2.0.CO;2

Holtslag, A. A. M., Svensson, G., Baas, P., Basu, S., Beare, B., Beljaars, A. C. M., et al. (2013). Stable boundary layers and diurnal cycles. *Bulletin America Meteorology Social, 94*(11), 1691–1706. https://doi.org/10.1175/BAMS-D-11-00187.1

Hourdin, F., Couvreaux, F., & Menut, L. (2002). Parameterization of the dry convective boundary layer based on a mass flux representation of thermals. *Journal of the Atmospheric Sciences, 59*(6), 1105–1123. https://doi.org/10.1175/1520-0469(2002)059<1105:potdcb>2.0.CO;2

Hourdin, F., Jam, A., Rio, C., Couvreux, F., Sandu, I., Lefebvre, M.-P., et al. (2019). Unified parameterization of convective boundary layer transport and clouds with the thermal plume model. *Journal of Advances in Modeling Earth Systems, 11*(9), 2910–2933. https://doi.org/10.1029/2019MS001666

Hourdin, F., Rio, C., Grandpeix, J.-Y., Madeleine, J.-B., Cheruy, F., Rochetin, N., et al. (2020). LMDZ6A: The atmospheric component of the IPSL climate model with improved and better tuned physics. *Journal of Advances in Modeling Earth Systems, 12*(7), e2019MS001892. https://doi.org/10.1029/2019MS001892

Hourdin, F., Williamson, D., Rio, C., Couvreux, F., Roehrig, R., Villefranque, N., et al. (2021). Process-based climate model development harnessing machine learning: II. model calibration from single column to global. *Journal of Advances in Modeling Earth Systems, 13*(6), e2020MS002225. https://doi.org/10.1029/2020MS002225

Jose Rodriguez-Manfredi. (2019). *APSS twins data*. NASA Planetary Data System. https://doi.org/10.17189/1518950
Kahre, M. A., Murphy, J. R., Newman, C. E., Wilson, R. J., Cantor, B. A., Lemmon, M. T., & Wolff, M. J. (2017). The Mars dust cycle. In *The atmosphere and climate of Mars* (pp. 295–337). Cambridge University Press. https://doi.org/10.1017/9781139060172.010

Knalay, E., & Kanamitsu, M. (1988). Time schemes for strongly nonlinear damping equations. *Monthly Weather Review*, 116(10), 1945–1958. https://doi.org/10.1175/1520-0493(1988)116(1945:TSFSND)2.0:CO;2

Kolmogorov, A. N. (1941). Energy dissipation in locally isotropic turbulence. *Doklady AN SSSR*, 32, 19–21.

Lange, L., Forget, F., Dupont, E., Vandemeulebrouck, R., Spiga, A., Millour, E., et al. (2023). Modeling slope microclimates in the Mars planetary climate model. *Journal of Geophysical Research: Planets*, 128(10), e2023JE007915. https://doi.org/10.1029/2023JE007915

Lenderink, G., & Holtslag, A. A. M. (2000). Evaluation of the kinetic energy approach for modeling turbulent fluxesin stratocumulus. *Monthly Weather Review*, 128(1), 244–258. https://doi.org/10.1175/1520-0493(2000)128(0244:EOTKEA)2.0:CO;2

Lenderink, G., & Holtslag, A. A. M. (2004). An updated length-scale formulation for turbulent mixing in clear and cloudy boundary layers. *Quarterly Journal of the Royal Meteorological Society*, 130(604), 3405–3427. https://doi.org/10.1256/qj.03.117

Li, D. (2019). Turbulent Prandtl number in the atmospheric boundary layer—Where are we now? *Atmospheric Research*, 216, 86–105. https://doi.org/10.1016/j.atmosres.2018.09.015

Li, D., Katul, G. G., & Zilitinkevich, S. S. (2016). Closure schemes for stably stratified atmospheric flows without turbulence cutoff. *Journal of the Atmospheric Sciences*, 73(12), 4817–4832. https://doi.org/10.1175/JAS-D-16-0101.1

Lopez-Gomez, I., Cohen, Y., He, J., Jaruga, A., & Schneider, T. (2020). A generalized mixing length closure for eddy-diffusivity mass-flux schemes of turbulence and convection. *Journal of Advances in Modeling Earth Systems*, 12(11), e2020MS002161. https://doi.org/10.1029/2020MS002161

Louis, J.-F. (1979). A parametric model of vertical eddy fluxes in the atmosphere. *Boundary-Layer Meteorology*, 17(2), 187–202. https://doi.org/10.1007/BF00712379

Louis, J. F., Tiedtke, M., & Geleyn, J.-F. (1982). A short history of the operational pbl parametrization at ECMWF. In *Paper presented at the ECMWF workshop on boundary layer parametrization, ECMWF, Reading*.

Marquet, P. (2008). Documentation of the TKE-CBR turbulent scheme in the CLIMATE version 5.0 of ARPEGE (based on the cycle 32t0 op1v2 13 + CNRM/GMGEC modset). In *GMAP documentation*. Retrieved from https://www.umr-cnrm.fr/gmapdoc/IMG/pdf_TKE_CBR_V2_eng.pdf

Mašek, J., Ďurán, I. B., & Brožková, R. (2022). Stable numerical implementation of a turbulence scheme with two prognostic turbulence energies. *Monthly Weather Review*, 150(7), 1667–1688. https://doi.org/10.1175/MWR-D-21-0172.1

Mellor, G. L., & Yamada, T. (1974). A hierarchy of turbulence closure models for planetary boundary layers. *Journal of the Atmospheric Sciences*, 30(6), 1061–1069. https://doi.org/10.1175/1520-0469(1973)030<1061:apotpo>2.0:co;2

Mellor, G. L., & Yamada, T. (1982). Development of a turbulence closure model for geophysical fluid problems. *Reviews of Geophysics and Space Physics*, 20(4), 851–875. https://doi.org/10.1029/rg020i004p00851

Montabone, L., Forget, F., Millour, E., Wilson, R., Lewis, S., Cantor, B., et al. (2015). Eight-year climatology of dust optical depth on Mars. *Icarus*, 251, 65–95. https://doi.org/10.1016/j.icanus.2014.12.034

Nieuwtsadt, F. T. M. (1984). The turbulent structure of the stable, nocturnal boundary layer. *Journal of the Atmospheric Sciences*, 41, 2202–2217.

Petrosyan, A., Galperin, B., Larsen, S. E., Lewis, S. R., Määttänen, A., Read, P. L., et al. (2011). The martian atmospheric boundary layer. *Reviews of Geophysics*, 49(3), RG3005. https://doi.org/10.1029/2010RG000351

Pottier, A., Forget, F., Montmessin, F., Navarro, T., Spiga, A., Millour, E., et al. (2017). Unraveling the Martian water cycle with high-resolution global climate simulations. *Icarus*, 291, 82–106. https://doi.org/10.1016/j.icanus.2017.02.016

Redelsperger, J.-L., Mahé, F., & Carlotti, P. (2001). A simple and general subgrid model suitable both for surface layer and free-stream turbulence. *Boundary-Layer Meteorology*, 101(3), 375–408. https://doi.org/10.1023/A:1019206001922

Rio, C., Hourdin, F., Couvreux, F., & Jam, A. (2010). Resolved versus parametrized boundary-layer plumes. Part II: Continuous formulations of mixing rates for mass-flux schemes. *Boundary-Layer Meteorology*, 135(3), 469–483. https://doi.org/10.1007/s10546-010-9478

Rodier, Q., Masson, V., Couvreux, F., & Paci, A. (2017). Evaluation of a buoyancy and shear based mixing length for a turbulence scheme. *Frontiers in Earth Science*, 5, 65. https://doi.org/10.3389/feart.2017.00065

Sandu, I., Beljaars, A., Bechtold, P., Mauritsen, T., & Balsamo, G. (2013). Why is it so difficult to represent stably stratified conditions in numerical weather prediction (NWP) models. *Journal of Advances in Modeling Earth Systems*, 5(2), 117–133. https://doi.org/10.1002/jame.20013

Spiga, A., Forget, F., Lewis, S. R., & Hinson, D. P. (2010). Structure and dynamics of the convective boundary layer on Mars as inferred from large-eddy simulations and remote-sensing measurements. *Quarterly Journal of the Royal Meteorological Society*, 136(647), 414–428. https://doi.org/10.1002/qj.563

Stull, R. B. (1990). *An introduction to boundary layer meteorology*. Kluver.

Sun, J. (2011). Vertical variations of mixing lengths under neutral and stable conditions during cases-99. *Journal of Applied Meteorology and Climatology*, 50(10), 2030–2041. https://doi.org/10.1175/JAMC-D-10-05006.1

Svensson, G., Holtslag, A. A. M., Kumar, V., Mauritsen, T., Steeneveld, G. J., Angevine, W. M., et al. (2011). Evaluation of the diurnal cycle in the atmospheric boundary layer over land as represented by a variety of single-column models: The second GABLS experiment. *Boundary-Layer Meteorology*, 140(2), 177–206. https://doi.org/10.1007/s10546-011-9611-7

Tiedtke, M. (1989). A comprehensive mass flux scheme for cumulus parameterization in large-scale models. *Monthly Weather Review*, 117(8), 1779–1800. https://doi.org/10.1175/1520-0493(1989)117(1779:ACMFSF)2.0:CO;2

van der Linden, S. J., Edwards, J. M., van Heerwaarden, C. C., Vignon, E., Genthon, C., Petenko, I., et al. (2019). Large-eddy simulations of the steady wintertime Antarctic boundary layer. *Boundary-Layer Meteorology*, 173(2), 165–192. https://doi.org/10.1007/s10546-019-00461-4

van de Wiel, B. J. H., Moene, A. F., Steeneveld, G. J., Baas, P., Bosveld, F. C., & Holtslag, A. A. M. (2010). A conceptual view on inertial oscillations and nocturnal low-level jets. *Journal of the Atmospheric Sciences*, 67, 2679–2689. https://doi.org/10.1175/2010JAS3289.1

van de Wiel, B. J. H., Moene, A. F. H. D. R. W., & Jonker, H. J. J. (2008). Local similarity in the stable boundary layer and mixing-length approaches: Consistency of concepts. *Boundary-Layer Meteorol*, 128, 103–116.

van de Wiel, B. J. H., Vignon, E., Baas, P., van Hooijdonk, I. G. S., van der Linden, S. J. A., van Hooft, J. A., et al. (2017). Regime transition in near-surface temperature inversions: A conceptual model. *Journal of the Atmospheric Sciences*, 74, 1057–1073. https://doi.org/10.1175/JAS-D-16-0180.1

Venayagamoorthy, S. K., & Stretch, D. D. (2010). On the turbulent Prandtl number in homogeneous stably stratified turbulence. *Journal of Fluid Mechanics*, 644, 359–369. https://doi.org/10.1017/S002211200999293X

Vignon, E., Hourdin, F., Genthon, C., Gallée, H., Bazile, E., Lefebvre, M.-P., et al. (2017). Parametrization of the boundary layer over the Antarctic plateau in a general circulation model: ID simulations against summertime observations at dome C. *Journal of Geophysical Research: Atmospheres*, 122(13), 6818–6843. https://doi.org/10.1002/2017JD026802
Vignon, E., Hourdin, F., Genthon, C., Van de Wiel, B. J. H., Gallée, H., Madeleine, J., & Beaumet, J. (2018). Modeling the dynamics of the atmospheric boundary layer over the Antarctic plateau with a general circulation model. *Journal of Advances in Modeling Earth Systems*, *10*(1), 98–125. https://doi.org/10.1002/2017MS001184

Vignon, E., Raillard, L., Genthon, C., Del Guasta, M., Heymsfield, A. J., Madeleine, J.-B., & Berne, A. (2022). Ice fog observed at cirrus temperatures at Dome C Antarctic plateau. *Atmospheric Chemistry and Physics*, *22*(19), 12857–12872. https://doi.org/10.5194/acp-22-12857-2022

Vignon, E., van de Wiel, B. J. H., van Hooijdonk, I. G. S., Genthon, C., van der Linden, S. J. A., van Hooft, J. A., et al. (2017). Stable boundary layer regimes at dome C, Antarctica: Observation and analysis. *Quarterly Journal of the Royal Meteorological Society*, *143*(704), 1241–1253. https://doi.org/10.1002/qj.2998

Weiss, B. P., & Ingersoll, A. P. (2000). Cold spots in the martian polar regions: Evidence of carbon dioxide depletion? *Icarus*, *144*(2), 432–435. https://doi.org/10.1006/icar.1999.6302

Yamada, T. (1983). Simulations of nocturnal drainage flows by a $q^2l$ turbulence closure model. *Journal of the Atmospheric Sciences*, *40*(1), 91–106. https://doi.org/10.1175/1520-0469(1983)040<0091:sondfb>2.0.co;2

Yamada, T., & Mellor, G. (1975). A simulation of the wangara atmospheric boundary layer data. *Journal of the Atmospheric Sciences*, *32*(12), 2309–2320. https://doi.org/10.1175/1520-0469(1975)032<2309:asotwa>2.0.co;2

Zilitinkevich, S., Elperin, T., Kleeorin, N., & Rogachevskii, I. (2007). Energy and flux-budget turbulence closure model for stably stratified flows. Part I: Steady state, homogeneous regime. *Boundary-Layer Meteorology*, *125*(2), 167–191. https://doi.org/10.1007/s10546-007-9189-2

Zilitinkevich, S., Elperin, T., Kleeorin, N., Rogachevskii, I., Esau, I., Mauritsen, T., & Miles, M. W. (2008). Turbulence energetics in stably stratified geophysical flows: Strong and weak mixing regimes. *Quarterly Journal of the Royal Meteorological Society*, *134*(633), 793–799. https://doi.org/10.1002/qj.264