# A Framework for Assessing Ocean Mixed Layer Depth Evolution

Alexandre Legay, Bruno Deremble, Thierry Penduff, Pierre Brasseur,

Jean-Marc Molines

► To cite this version:

Alexandre Legay, Bruno Deremble, Thierry Penduff, Pierre Brasseur, Jean-Marc Molines. A Framework for Assessing Ocean Mixed Layer Depth Evolution. Journal of Advances in Modeling Earth Systems, 2024, 16 (10), pp.e2023MS004198. 10.1029/2023MS004198. hal-04727426

HAL Id: hal-04727426

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

Submitted on 9 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2023MS004198

**Key Points:**

* A parameter space is proposed to assess the evolution of the mixed layer depth for realistic forcings and preconditioning conditions

* An evaluation of a collection of 1D simulations shows a statistically good performance of the parameter space

* Two applications demonstrate the utility of the parameter space for assessing and comparing realistic 3D simulations

**Correspondence to:**

A. Legay,
alexandre.legay@univ-grenoble-alpes.fr

**Citation:**

Legay, A., Deremble, B., Penduff, T., Brasseur, P., & Molines, J.-M. (2024). A framework for assessing ocean mixed layer depth evolution. *Journal of Advances in Modeling Earth Systems, 16*, e2023MS004198. https://doi.org/10.1029/2023MS004198

Received 22 DEC 2023

Accepted 24 JUL 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_2.png" /></figure>

# A Framework for Assessing Ocean Mixed Layer Depth Evolution

Alexandre Legay¹ ID, Bruno Deremble¹ ID, Thierry Penduff¹ ID, Pierre Brasseur¹, and Jean-Marc Molines¹

¹CNRS, INRAE, IRD, Grenoble INP, IGE, Université Grenoble Alpes, Grenoble, France

**Abstract** The ocean surface mixed layer plays a crucial role as an entry or exit point for heat, salt, momentum, and nutrients from the surface to the deep ocean. In this study, we introduce a framework to assess the evolution of the mixed layer depth (MLD) for realistic forcings and preconditioning conditions. Our approach involves a physically-based parameter space defined by three dimensionless numbers: $\lambda_s$ representing the relative contribution of the buoyancy flux and the wind stress at the air-sea interface, $R_h$ the Richardson number which characterizes the stability of the water column relative to the wind shear, and $f/N_h$ which characterizes the importance of the Earth's rotation (ratio of the Coriolis frequency $f$ and the pycnocline stratification $N_h$). Four MLD evolution regimes ("restratification," "stable," "deepening," and "strong deepening") are defined based on the values of the normalized temporal evolution of the MLD. We evaluate the 3D parameter space in the context of 1D simulations and we find that considering only the two dimensions ($\lambda_s$, $R_h$) is the best choice of 2D projection of this 3D parameter space. We then demonstrate the utility of this two-dimensional $\lambda_s - R_h$ parameter space to compare 3D realistic ocean simulations: we discuss the impact of the horizontal resolution (1°, 1/12°, or 1/60°) and the Gent-McWilliams parameterization on MLD evolution regimes. Finally, a proof of concept of using observational data as a truth indicates how the parameter space could be used for model calibration.

**Plain Language Summary** Vertical mixing of water near the ocean surface occurs when cold air temperatures create dense cold water at the surface that tends to sink in the ocean or when a strong wind induces turbulence at the ocean surface. These processes mix heat and salt and create a layer at the top of the ocean that has a uniform temperature and salinity and that is called the “mixed layer.” This mixed layer plays a fundamental role in the Earth climate system, and the representation of its evolution in ocean models hence needs to be assessed. For this purpose, we propose to map the mixed layer evolution in a three-dimensional space where the first axis is related to the wind and the surface heat flux, the second axis to the stability of the water column, and the third axis to the Earth's rotation. We show that this tool performs statistically well and we present how to use it in the context of realistic ocean models.

## 1. Introduction

The evolution of the mixed layer near the air-sea interface is primarily driven by the vertical mixing and restratification processes. Vertical mixing is usually driven by winds, surface cooling, brine rejection, Langmuir turbulence, and wave breaking (Q. Li et al., 2019; Marshall & Schott, 1999; Vreugdenhil & Gayen, 2021). In contrast, restratification processes are driven by solar heating, freshwater flux, or lateral processes such as mixed layer instabilities (see for example Boccaletti et al., 2007; Fox-Kemper et al., 2007). Accurately representing the mixed layer depth (MLD) evolution is crucial for capturing many physical and biogeochemical mechanisms, such as the sequestration of heat and carbon by the ocean (e.g., Banks & Gregory, 2006; Bernardello et al., 2014), the dynamics of marine ecosystems (e.g., Lévy et al., 1998; Sverdrup, 1953; Taylor & Ferrari, 2011), and the representation of the Atlantic meridional overturning circulation (e.g., Kuhlbrodt et al., 2007).

Historically, various approaches have been proposed to describe the vertical mixing in the MLD which can be grouped into three main categories: bulk mixed layer models, similarity models, and turbulence closure models. For bulk mixed layer models, the governing equations of fluid dynamics are integrated over the mixed layer and represent the evolution of integrated properties (e.g., Gaspar, 1988; Kraus & Turner, 1967; Pollard et al., 1973; Price et al., 1986). These models have been used to derive theoretical scalings for the evolution of the MLD, such as the wind-driven deepening $h \propto u_* N^{-1/2} t^{1/2}$ (Pollard et al., 1973), observed empirically by Price (1979), and the
free convection scaling $h \propto Q^{1/2} N^{-1} t^{1/2}$ (Turner, 1973; Van Roekel et al., 2018) measured empirically by Souza et al. (2020) ($h$ being the MLD, $u_*$ the surface friction velocity, $t$ the time, $Q$ the net surface heat flux and $N$ the Brunt Väisälä frequency).

The second class of models are the similarity models such as the K-Profile Parameterization (KPP, Large et al., 1994) or the OSMOSIS model (Damerell et al., 2020; Madec et al., 2022). These models assume that the vertical profiles of tracers and momentum are self-similar. With this self-similarity hypothesis, turbulent fluxes can be computed by scaling a predefined profile shape with the magnitude of the surface forcing. Although the KPP model successfully captures many observed features of the ocean's boundary layer, it relies on empirical relationships and is not derived from first principles. Nevertheless, KPP remains one of the most widely used parameterizations of vertical mixing in ocean models (Q. Li et al., 2019; Souza et al., 2020; Van Roekel et al., 2018).

The last class of models consists of the turbulence closure models. These models consider equations of higher order moments of the turbulent quantities and make some assumptions about their formulations in order to close the problem (e.g., Mellor, 1973; Mellor & Yamada, 1974, 1982). Widely used models in this class include the Turbulent Kinetic Energy (TKE) models, which solve a prognostic equation for TKE (Gaspar et al., 1990), and the Generic Length Scale (GLS) models which include an additional prognostic equation for a turbulence length scale quantity (global description: Umlauf & Burchard, 2003, 2005; examples of models of this type: $k-\epsilon$: Hanjalić & Launder, 1972; Rodi, 1987, $k-kl$: Mellor & Yamada, 1982, $k-\omega$: Wilcox, 1988, $k-\tau$: Zeierman & Wolfstein, 1986; Thangam et al., 1992).

Some models do not fit into one of the three aforementioned classes, such as the energetics-based Planetary Boundary Layer scheme (ePBL, Reichl & Hallberg, 2018), which combines a depth-dependent bulk mixed layer model with a turbulence closure model.

Currently, most climate simulations either use TKE or KPP models for vertical mixing (Zhu et al., 2020). The MLD evolution of these climate simulations depends on (a) the choice of the vertical mixing scheme, (b) the impact of resolved lateral processes, and (c) the use of parameterizations for unresolved lateral effects (e.g., the Fox-Kemper et al. (2007) and the Gent and McWilliams (1990) parameterizations). Consequently, objectively comparing the MLD evolution in these climate simulations is challenging (Treguier et al., 2023). Common approaches involve comparing hydrographic sections at specific locations (e.g., evaluation at the Papa station: Burchard & Bolding, 2001; Gaspar et al., 1990; Giordani et al., 2020; Large et al., 1994), and/or conducting intercomparisons at specific times (e.g., with intercomparisons of MLD maps: Gutjahr et al., 2021; Heuzé, 2017), and/or using indirect metrics of the MLD evolution (e.g., by comparing the amount of deep water formed: Koenig et al., 2021). However, these approaches only explore a limited range of forcings and preconditioning conditions.

In this paper, we adopt a more comprehensive approach by considering all possible ranges of forcings and preconditioning conditions in a suitable parameter space. Belcher et al. (2012) and Q. Li et al. (2019) have pioneered this approach to evaluate the representation of Langmuir circulation with different vertical mixing schemes. They proposed a two-dimensional parameter space, with the first dimension (i.e., the first dimensionless number) assessing the relative importance of the wind and the wave forcings, and the second dimension characterizing the relative importance of wave and buoyancy forcings. In this parameter space, Belcher et al. (2012) and Q. Li et al. (2019) defined theoretical boundaries to highlight the importance of different surface forcings. Subsequently, Large eddy simulations (LES) results found in the literature were placed in this parameter space to see which regimes are explored by these simulations. Their objective was to identify a potential bias arising from miscalibration of the LES simulations used to establish parameterizations.

Following the approach of Q. Li et al. (2019), our objective is to propose a dedicated parameter space to describe the evolution of the MLD. This parameter space aims to capture the MLD evolution dependency on the relative importance of wind and buoyancy forcings, preconditioning conditions (Marshall & Schott, 1999), and the influence of Earth's rotation. More precisely, we will evaluate the relative deepening or shoaling of the MLD over a 1-day period ($\partial_t h/h$ from noon to noon expressed in %/day). To keep the practicability of having few parameters, we have decided to exclude several processes (such as waves). Our study demonstrates that, at first order, three dimensionless numbers are enough for characterizing MLD evolution. In contrast to Q. Li et al. (2019), our
approach involves directly plotting the values of MLD evolution $\partial_{h}h/h$ in the parameter space. This direct visualization allows for a more straightforward comparison of the behavior of different simulations.

This article is constructed as follows. First, we present in Section 2 the three dimensionless numbers that constitute the parameter space. Second, we show in Section 3.1 that MLD evolution regimes naturally emerge in this parameter space in the context of 1D simulations. Third, we present in Sections 3.2 and 3.3 two applications for showing how the parameter space can be used in practice with 3D realistic ocean models. The first application is about the impact of the lateral resolution on the MLD evolution regimes. The second one focuses on the effect of the Gent McWilliams (GM) parameterization which aims at representing the impact of the unresolved mesoscale processes in a coarse-resolution ocean model. Finally, we conclude in Section 4 on the practical use of this three-dimensional parameter space, discuss its strengths and limitations, and open a door about using observations as a truth in the parameter space to aid model calibrations.

## 2. Materials and Methods

### 2.1. Definition of the Three Dimensionless Numbers

We will formulate three dimensionless numbers to characterize the evolution of the MLD $h$ of a 1D water column model that evolves from one day to another according to daily-mean surface forcings and preconditioning conditions. The goal here is to identify the main factors that drive the MLD evolution $\partial_{h}h$ in order to build a parameter space with a small number of dimensions being useable for evaluating $\partial_{h}h$. The principal omissions resulting from our choices will be discussed further in Section 2.2.

We adopt the description of the water column given by bulk mixed layer models (e.g., Pollard et al., 1973): near the surface, we consider a well-mixed layer of thickness $h$. This layer is forced at the surface by the wind stress with a friction velocity $u_{*} = [\overline{(u'w')_{z=0}}]^2 + (\overline{v'w'}_{z=0})^2]^{1/4}$ and the downward surface buoyancy flux $B_0 = -\overline{w'b'}_{z=0}$ ($B_0 < 0$ for a destabilizing flux at the ocean surface), with $u'$, $v'$, and $w'$ the turbulent velocities, the overline that denotes an average over small scale fluctuations (see Stull, 1988), $b'$ the fluctuation of the buoyancy $b = -\frac{\rho-\rho_0}{\rho_0}g$, $\rho$ the density, $\rho_0$ the reference density, and $g$ the acceleration due to gravity. At the base of the mixed layer, the stratification is given by the Brunt Vaisala frequency $N_h$, which is sometimes called “pre-conditioning.” In order to describe the MLD evolution, we have also opted to retain the local Coriolis parameter $f$. With this idealized view of the mixed layer, the MLD evolution $\partial_{h}h$ is a function of five physical quantities: ($u_*$, $B_0$, $h$, $N_h$, $f$). These five physical quantities are expressed with 2 distinct dimensions: length and time ($[u_*] = L T^{-1}$, $[B_0] = L^2 T^{-3}$, $[h] = L$, $[N_h] = T^{-1}$, $[f] = T^{-1}$). The Vaschy-Buckingham theorem ($\pi$ theorem) thus states that these five physical quantities can be represented by $5 - 2 = 3$ dimensionless numbers. The three dimensionless numbers we have chosen are

$$ \lambda_s = \frac{-B_0 h}{u_*^3}, \qquad (1) $$

$$ R_h = \left( \frac{N_h h}{u_*} \right)^2, \qquad (2) $$

and

$$ f/N_h. \qquad (3) $$

Note that $\lambda_s$ is positive for a destabilizing surface buoyancy flux ($B_0 < 0$) and negative for a stabilizing surface buoyancy flux ($B_0 > 0$).

Henceforth, we describe the physical interpretations of these three dimensionless numbers and then present the associated three-dimensional parameter space.

#### 2.1.1. Physical Interpretation of $\lambda_s$

In the context of ocean mixed layer dynamics, $\lambda_s$ can be interpreted in at least three ways.
*First Interpretation:* $\lambda_s$ can be interpreted using the evolution equation of the TKE (for a full description of this equation, see Hanjalić & Launder, 1972; Rodi, 1987; Umlauf & Burchard, 2003):

$$ \frac{Dk}{Dt} = P + G - \epsilon + D_k \quad (4) $$

with $k = \frac{1}{2}(\overline{u^2} + \overline{v^2} + \overline{w^2})$ the TKE, $P = -\overline{u'w'} \partial_z \overline{u} - \overline{v'w'} \partial_z \overline{v}$ the TKE production by the shear (by conversion of mean kinetic energy), $G = \overline{w'b'}$ the TKE production (or destruction) by the turbulent buoyancy flux (by conversion of mean potential energy), $\epsilon$ the TKE dissipation, and $D_k$ the TKE diffusion.

The surface layer is defined as the zone where the turbulent fluxes vary by less than 10% from their values at $z = 0$ (Stull, 1988). In this zone, we can consider $\overline{u'w'} \approx \overline{u'w'}|_{z=0} \propto u_*^2$ and $G \approx \overline{w'b'}|_{z=0} = -B_0$. Since the surface mean horizontal velocity $\overline{u}|_{z=0}$ is well correlated to the surface friction velocity $u_*$ (Weber, 1983), and if we neglect the mean horizontal velocity $\overline{u}$ below the MLD (Pollard et al., 1973), then an order of magnitude of $\partial_z \overline{u}$ is given by $u_*/h$. It follows $P \propto u_*^3/h$. An evaluation of $G/P$ in the surface layer finally gives

$$ \left. \frac{G}{P} \right|_{surf} \propto \frac{-B_0}{\left(\frac{u_*^3}{h}\right)} \equiv \lambda_s. \quad (5) $$

This ratio $G/P|_{surf}$ is by definition the flux Richardson number $R_f$ evaluated in the surface layer (Mellor & Durbin, 1975). It gives the relative contribution of surface buoyancy flux and wind for the production of TKE. In the case $B_0 < 0$ that is, $G/P|_{surf} > 0$, both terms produce TKE. On the other hand, in the case of a restratifying buoyancy flux $B_0 > 0$, there is a competition between production by the shear ($P > 0$) and destruction by the turbulent buoyancy flux ($G < 0$). Particularly, for $G/P|_{surf} < -1$, more TKE is destroyed (converted into mean potential energy) than created (from mean kinetic energy): this likely represents a restratification event.

*Second Interpretation:* We can also interpret $\lambda_s$ in the light of the Monin-Obukhov similarity theory (Obukhov, 1971). This theory, which is valid in the surface layer, introduces the Monin-Obukhov length $L_{MO}$:

$$ L_{MO} = \frac{u_*^3}{\kappa B_0}. \quad (6) $$

We give here its definition in the oceanic framework (see for example Zheng et al., 2021) which is the opposite of the atmospheric definition. The physical interpretation of $L_{MO}$ was introduced by Obukhov in the case $L_{MO} < 0$ ($\Leftrightarrow B_0 < 0 \Leftrightarrow \lambda_s > 0$). In this regime, $L_{MO}$ estimates the typical thickness of a “sub-layer of dynamic turbulence” in which stratification is of little importance and the turbulence dynamics is governed by the mean-current shear (Obukhov, 1971), that is, the production of TKE by the buoyancy $G$ is negligible in comparison to the one by the mean-current shear $P$. In practice, Wyngaard (1973) has shown that $G \approx P$ for $z \approx 0.5L_{MO}$ (Figure 5.22, Stull, 1988). The number $\lambda_s$ can be seen as

$$ \lambda_s = \frac{1}{\kappa} \frac{-h}{L_{MO}}. \quad (7) $$

Thus, $\lambda_s < 0.5/\kappa$ gives $P > G$ in the mixed layer while $\lambda_s > 0.5/\kappa$ means $G > P$. It is important to recall that this interpretation only stands for $L_{MO} < 0$ ($\Leftrightarrow B_0 < 0 \Leftrightarrow \lambda_s > 0$). For the case $L_{MO} > 0$ ($\Leftrightarrow B_0 > 0 \Leftrightarrow \lambda_s < 0$), we refer to our first interpretation of $\lambda_s$.

*Third Interpretation:* In the case $B_0 < 0$, convective thermals have a velocity in the order of $w_* = (-B_0h)^{1/3}$ (Marshall & Schott, 1999; Willis & Deardorff, 1974). Then $\lambda_s$ can be written as

$$ \lambda_s = \left( \frac{w_*}{u_*} \right)^3. \quad (8) $$

In this expression, it is clear that $\lambda_s$ measures the relative importance of mechanical and convective forcings.
Last, it is worth noting that in a different context Simpson and Hunter (1974) used a similar ratio to characterize the mixing occurring in the Irish Sea where $u_*$ was related to the tidal forcing (friction in the bottom boundary layer).

### 2.1.2. Physical Interpretation of $R_h$

The dimensionless number $R_h$ can be interpreted as a Richardson number. By definition the gradient Richardson number $R_i = N^2/(\partial_z \bar{u})^2$ is the ratio of the stabilizing effect of the stratification and the destabilizing effect of the shear of the mean current (see for example Mack & Schoeberlein, 2004). We compare here the stratification at the mixed layer base $N_h^2$ with the order of magnitude of the wind-induced shear $u_*/h$. This gives the ratio:

$$ \frac{N_h^2}{\left(\frac{u_*}{h}\right)^2} \equiv R_h. \qquad (9) $$

We could have included the contribution of $w_*$ to the definition of the shear but we will see in Section 2.3 that this omission is intentional and results in a simpler interpretation of the parameter space.

### 2.1.3. Physical Interpretation of $f/N_h$

In this paragraph, we endeavor to keep the discussion fairly general in scope, even if the numerical models that we use in this study do not resolve all the processes described as we shall see later. There are several interpretations for the dimensionless number $f/N_h$ in the context of mixed layer dynamics.

A well known effect of rotation is that wind triggers inertial oscillations which slow down the mixed layer deepening when $h > h_p$, with $h_p \propto u_*/\sqrt{N_h f}$ the maximum MLD predicted by Pollard et al. (1973) in case of a shear-driven mixed layer deepening in a rotating case. To isolate the effect of rotation and highlight the non-dimensional number $f/N_h$, we can write

$$ \frac{h}{h_p} = \sqrt{R_h} \sqrt{\frac{f}{N_h}}. \qquad (10) $$

At finer scales, rotation also plays a key role in the organization of submesoscale dynamics and mixed layer instabilities (Taylor & Thompson, 2023). The ratio $f/N_h$ can be seen as $h/L_d$, where $L_d \propto N_h h/f$ is the “mixed layer” Rossby radius of deformation in the quasi-geostrophic context ($h$ is not the total depth of the fluid but the MLD). In the situation where $f/N_h > O(1)$, we expect that mixed layer instabilities will create a lateral buoyancy flux (see Boccaletti et al., 2007). Part of the turbulent energy normally used for vertical mixing is hence used for lateral mixing. Therefore, we expect that values $f/N_h > O(1)$ result in a slowdown of the MLD deepening.

At even finer scales, rotation drives the formation of small vortices in convective plumes which slow down the MLD deepening (Deremble, 2016; Dingwall & Taylor, 2024; Frank et al., 2021), so we expect that for high values of $f/N_h$, we are likely to observe a reduced MLD deepening (a phenomenon that is reminiscent of the rapidly rotating columnar regime for which we observe a decrease of the turbulent heat flux; see Bouillaut et al., 2019; Aurnou et al., 2020). Speer and Marshall (1995) have also described how the aspect ratio of convective plumes is determined by the ratio $f/N_h$, where the effect of rotation is mainly to alter the lateral spreading of convective structures.

## 2.2. Limitations

There are of course other physical phenomena that occur in the mixed layer and that we have not taken into account:

* The effect of waves and associated Langmuir turbulence that could have been represented through the values of the surface Stokes drift $u_0^S$ (Q. Li et al., 2019). However, it is worth mentioning that part of $u_0^S$ can be explained by $u_*$. Minimal parameterizations of Langmuir turbulence even define $u_0^S$ directly proportional to $u_*$ (M. Li & Garrett, 1993; Madec et al., 2022). Thus, some of the wave impacts are implicitly contained through the consideration of $u_*$.
* All the effects of the horizontal gradients (of velocities, pressure...) and advections that are present in a 3D realistic ocean model. We can particularly pinpoint the Ekman flow that can create an equivalent stabilizing/ destabilizing wind-driven buoyancy flux (see for example Thomas & Lee, 2005), and the impacts of the restratification by baroclinic instability at convective fronts or at mesoscale eddy fronts. Some of these aspects can be captured by looking at the isopycnal slopes and this point will be further discussed in Section 3.3.

* The influence of the vertical shear of the horizontal velocities. As we work with daily evolution, and knowing the Ekman theory (Ekman, 1905), we can expect this shear to be partly represented by the consideration of the surface wind friction velocity, the MLD, and the Coriolis parameter.

Considering many of these aspects would have meant adding more dimensions to the parameter space and thus reducing its practical use. We will see in the results that the dimensionless numbers we have chosen are in many situations sufficient to obtain a significant prediction of the MLD evolution and therefore capture well the dominant processes of this evolution.

## 2.3. Other Possible Dimensionless Numbers

We have defined three dimensionless numbers from the five physical quantities ($u_*$, $B_0$, $h$, $N_h$, $f$), but other choices could have been possible. We want to highlight here some of them:

* The Rossby number $Ro = u_*/(hf)$ that characterizes the relative importance of the inertial and the Coriolis forces (Van Der Laan et al., 2020).

* The ratio $h/h_p$ compares the depth $h$ to the maximum one $h_p$ predicted by Pollard et al. (1973). In idealized cases with only wind and no buoyancy flux, this dimensionless number is a good predictor of the evolution of the mixed layer (Ushijima & Yoshikawa, 2020). However, in realistic conditions, we found that this dimensionless number does not add much information compared to $R_h$ (not shown).

* In case of $B_0 < 0$, the Richardson number $R_h^* = (N_h h/w_*)^2$ constructed with $w_*$ rather than $u_*$ which compares the stabilizing effect of the stratification and the destabilizing impact of the buoyancy flux (Shy, 1995; Turner, 1986). However, $R_h^*$ can be expressed as

$$R_h^* = R_h/\lambda_s^{2/3} \qquad (11)$$

and hence, in a log-log parameter space ($\lambda_s$, $R_h$), we will see that the isolines of $R_h^*$ appear as lines of slope 2/3. Then, we have the possibility to see the isolines of both $R_h$ and $R_h^*$ at the same time in the parameter space. We will use this information to decide whether one or the other is more representative of the deepening situation. Such distinction would have been impossible if we had taken a Richardson number defined with the two contributions at the same time, such as $(N_h h/\sqrt{u_* w_*})^2$, $(N_h h/\max(u_*, w_*))^2$ or $(N_h h/(u_* + w_*))^2$.

The next four sections present the simulations we will use to evaluate the parameter space and to conduct the two model sensitivity studies. In the next three sections, we present three 3D realistic ocean simulations at three different horizontal resolutions (1/60°, 1/12°, and 1°). These simulations will be used to study the impact of the horizontal resolutions on the MLD evolution regimes. Table 1 summarizes the main features of these simulations.

## 2.4. NEMO-eNATL60 1/60° Basin-Scale North Atlantic Ocean Simulation

The eNATL60-BLBT02 (eNATL60) simulation (Brodeau et al., 2020) is a basin-scale North Atlantic ocean/sea-ice simulation forced by the atmospheric three-hourly ERA-Interim reanalysis (Dee et al., 2011) on a 1/60°-horizontal and 300-vertical-level grid (vertical resolution of 1 m near the surface and 50 m in the deep ocean). It includes an explicit tidal forcing. The lateral boundary conditions for the ocean velocities, temperature, and salinity are based on the GLORYS12 v1 reanalysis (Lellouche et al., 2021). Vertical mixing is governed by a TKE scheme combined with the Enhanced Vertical Diffusivity (EVD) parameterization which increases the vertical diffusivity in case of unstable water columns (Lazar et al., 1999; Madec et al., 2022). The Fox-Kemper parameterization (Fox-Kemper et al., 2007), which represents the restratifying effect of sub-mesoscale mixed layer eddies, is included. The simulation covers 10 months from 1 January 2010 to 29 October 2010. A description of the technical choices and the configuration files are available at https://github.com/ocean-next/eNATL60. For our study, we extracted daily averages of the data in two regions of interest (Figure 1). The “Western Mediterranean region” extends from 2°E to 10°E and from 40°N to 44°N, and the “Labrador region” from 56°W to 51°W
<table><caption>**Table 1**
Summary of the Main Features of the Simulations</caption><thead><tr><th rowspan="2">Simulation</th><th colspan="2">Horizontal resolution (subsampling)</th><th rowspan="2">Number of vertical levels (in 0 =-&gt;2,000 m)</th><th colspan="2">Total time</th><th rowspan="2">Vertical mixing scheme</th><th rowspan="2">Additional parameterization</th><th rowspan="2">Reference</th></tr><tr><th>Med/Lab</th><th>Global</th><th>Med/Lab</th><th>Global</th></tr></thead><tbody><tr><td>eNATL60</td><td>1/60° (1/12°)</td><td>—</td><td>300 (157)</td><td>10 months (January to October 2010)</td><td></td><td>TKE + EVD</td><td>Fox-Kemper</td><td>Brodeau et al. (2020) (https://github.com/ocean-next/eNATL60)</td></tr><tr><td>eORCA12</td><td>1/12°</td><td>(15°)</td><td>75 (54)</td><td>10 years (2006–2015)</td><td></td><td>TKE + EVD</td><td>Internal Wave Mixing (IWM)</td><td>https://github.com/meom-configurations/eORCA12.L75-GJM2020</td></tr><tr><td>eORCA1</td><td>1°</td><td>(15°)</td><td>75 (54)</td><td>20 years (2000–2019)</td><td></td><td>TKE + EVD</td><td>—</td><td>https://github.com/meom-configurations/eORCA1-GJM2020</td></tr><tr><td>eORCA1GM</td><td>1°</td><td>(15°)</td><td>75 (54)</td><td>20 years (2000–2019)</td><td></td><td>TKE + EVD</td><td>Gent McWilliams (GM)</td><td>https://github.com/meom-configurations/eORCA1-GJM2020</td></tr><tr><td>Collection of 1D simulations</td><td colspan="2">Simulations performed at the locations of the 1° grid</td><td>54</td><td>10 years (2006–2015)</td><td>20 years (2000–2019)</td><td>TKE + EVD</td><td>—</td><td>Legay (2024b) (https://github.com/legaya/James2024ParameterSpace/)</td></tr></tbody></table>

Note. “Med” and “Lab” stand for respectively the “Western Mediterranean” and the “Labrador” extractions (see text).

and from 55°N to 59°N. These regions were chosen because they are well known sites where moderate to intense convective events occur. To avoid shallow water coastal dynamics, we only kept locations for which the local depth is greater than or equal to 2,000 m. Moreover, to reduce the amount of data, we subsampled the horizontal resolution of the outputs from 1/60° to 1/12°. For doing that, we used the function “samplegrid” of the Climate Data Operators library (CDO; Schulzweida, 2023) with a subsampling factor of 5 on both x and y dimensions of the grid.

## 2.5. NEMO ORCA 1/12° Global Ocean Simulation

The eORCA12.L75-GJM2020 (eORCA12) simulation is a global ocean/sea-ice simulation forced by the atmospheric reanalysis JRA55-do 1.4.0 (Tsujino et al., 2018) performed over the period 1979–2019. It uses the ORCA12.L75 grid which has a 1/12° horizontal resolution and a 75-level non-uniform vertical grid (vertical resolution of 1 m near the surface and 200 m in the deep ocean). Vertical mixing is governed by a TKE + EVD + Internal Wave Mixing (additional parameterization accounting for mixing due to internal waves) scheme. A description of the technical choices and the configuration files are available at https://github.com/meom-configurations/eORCA12.L75-GJM2020. In this study, we used a 10-year period (1 January 2006 to 1 January 2015) and we extracted daily averages of the data on the same Western Mediterranean and Labrador regions as described above (Figure 1), with the same restriction of keeping only locations for which the local depth is greater than or equal to 2,000 m. Unlike eNATL60, which is not a global simulation, we are also going to use the results of eORCA12 at the global scale and, to reduce the amount of data, we subsampled the outputs from the 1/12° resolution to only keep 234 points placed on a 15° grid (cf. Figure 1). This coarse representation is enough to capture a realistic range of $f/N_h$, representative of the world's oceans.

## 2.6. NEMO ORCA 1° Global Ocean Simulation

We performed the eORCA1-GJM2020 (eORCA1) simulation in a set-up that is the same as the eORCA12 simulation, except for the horizontal resolution. However, a 1° horizontal resolution is not considered eddy-resolving and for this reason, an eddy-induced velocity is often added at this resolution to parameterize missing mesoscale eddies. We computed this eddy-induced velocity with the Gent McWilliams (GM) parameterization (Gent & McWilliams, 1990). We performed two experiments, eORCA1 without GM and eORCA1GM with GM. The GM coefficient was taken constant (with the NEMO default input parameters of lateral diffusive velocity $U_e = 0.02 \text{ m s}^{-1}$ and lateral diffusive length $L_e = 2 \times 10^5 \text{ m}$). A description of the technical choices and
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The three regions of interest in this study. (a, b) are respectively the “Western Mediterranean” region that extends from 2°E to 10°E and from 40°N to 44°N, and the “Labrador” region that extends from 56°W to 51°W and from 55°N to 59°N. (c) presents the 15° grid used at the global scale (Mollweide's projection). All of these three regions are restricted to depths greater than 2,000 m.</figcaption></figure>
the configuration files are available at https://github.com/meom-configurations/eORCA1-GJM2020. We took daily averages of the outputs over 20 years from 1 January 2000 to 1 January 2019. We extracted the outputs on the two regions of interest, as well as at the global scale subsampled on the 15° grid described in the previous section. Again, only the locations for which the local depth is greater than or equal to 2,000 m were kept.

## 2.7. 1D Simulations

This section presents the 1D simulations we will use to classify the MLD evolution regimes in the parameter space.

We performed a collection of 1D water column simulations using the code presented in Fearon et al. (2020), which is a standalone 1D vertical version of the Coastal and Regional Ocean COmmunity model (CROCO, https://www.croco-ocean.org/). This collection contains 20-year simulations from 1 January 2000 to 1 January 2019 at all the 1° grid locations of the two regions Labrador and Western Mediterranean, as well as 10-year simulation from 1 January 2006 to 1 January 2015 for all the points of the global subsampled grid of 15° resolution. To be consistent with the 3D models, we used a TKE + EVD scheme, we included the Earth's rotation, we kept only the locations where the local depth is greater than 2,000 m and we applied the same atmospherical forcings as the eORCA1 simulation (presented in Section 2.6). The wind, precipitation, evaporation, and non-solar heat flux forcings were applied daily with 24-hr constant values. The solar flux was constructed with a cosine truncated of its negative values, thus representing 12 hr daytime with positive values and 12 hr nighttime with zero values. For simplicity, the length of day is taken constant over the year (the seasonal cycle of the solar flux is however well maintained because it is indirectly contained in the values of the daily means). The temperature and salinity profiles were re-initialized to their eORCA1 values at the beginning of each new year, so these simulations should be viewed as multiple annual simulations. The vertical grid was taken equal to the one of eORCA1/eORCA12 cut at 2,000 m depth, hence the 54 shallowest levels of this 75-level grid. The time step was set to 360 s. The UNESCO 1983 nonlinear equation of state was used (Fofonoff & Millard, 1983). At the bottom boundary of the domain, we imposed a homogeneous Neumann condition (no flux).

## 2.8. Practical Calculation of the Dimensionless Numbers

The five physical quantities ($u_*$, $B_0$, $h$, $N_h$, $f$) appearing in the dimensionless numbers are calculated as follows:

* The surface wind friction velocity is calculated from its definition $u_* = \sqrt{|\tau^w|/\rho_0}$ with $|\tau^w|$ the norm of the wind stress vector at the ocean surface.

* We calculate the surface buoyancy flux with its classical linear definition: $B_0 = \frac{g}{\rho_0}(\frac{\alpha Q}{c_p} - \beta S_{surf}(E - P))$ with $Q$ the downward surface heat flux ($Q < 0$ for cooling), $E$ the evaporation, $P$ the precipitation, $c_p$ the heat capacity per unit mass, $S_{surf}$ the surface salinity, $\rho_0$ the reference density, $\alpha$ the thermal expansion coefficient, and $\beta$ the haline contraction coefficient. These two coefficients ($\alpha, \beta$) are calculated with the non-linear equation of state for each location at every time step according to the local values of the surface temperature and the surface salinity.

* We choose the MLD definition of the CMIP6 working group for $h$ (Griffies et al., 2016, Appendix H24.2). This definition is based on a buoyancy difference from the surface and was designed to give results similar to the density criterion of a 0.03 kg/m³ difference of de Boyer Montégut et al. (2004) in the case of a local density close to 1,035 kg/m³.

* The stratification at the base of the mixed layer $N_h^2$ is defined as a difference between the properties at the MLD ($z = -h$) and 10 % below the MLD ($z = -1.1$ hr):

$$N_h^2 = \frac{g \rho(z = -1.1h) - \rho(z = -h)}{0.1h} \qquad (12)$$

It is worth mentioning here that, if we interpret $\Delta\rho = \rho(z = -1.1 \text{ hr}) - \rho(z = -h)$ as the “density jump” at the base of the mixed layer, $R_h = (N_h h/u_*)^2$ can be written as $R_h = g \frac{\Delta\rho}{\rho_0} h/(0.1 u_*^2)$ and is hence proportional to $R_\tau = g \frac{\Delta\rho}{\rho_0} h/u_*^2$ the bulk Richardson number associated to the wind (Price, 1979).

* The Coriolis parameter is equal to $f = 2\Omega_0 \sin(\phi)$ with $\Omega_0 = 7.29 \times 10^{-5} \text{ rad s}^{-1}$ the rotation rate of the Earth and $\phi$ the latitude. As we did not consider any lateral gradients, the sign of $f$ should not matter for the scaling
and we thus consider the absolute value of f. However, throughout the manuscript, we write f rather than |f| for brevity.

Sensitivity of the results to other choices of MLD definitions (de Boyer Montégut et al., 2004; Reichl et al., 2022) and other N<sub>h</sub><sup>2</sup> evaluations (centered at z = -h or with a constant distance of 15 m below the mixed layer, see the discussion in Sérazin et al., 2023) were tested but not shown here for brevity. In short, the two definitions we chose were the ones giving the results with the highest significance (the notion of “significance” in the parameter space is defined in the next section) and hence the ones that are the more relevant in our context of the evaluation of the relative MLD deepening or shoaling over a 1-day period.

As we follow the MLD evolution over a 1-day period, we use the daily averages of the quantities (from noon to
noon). We opt to take the daily averages of *h* and *N*<sub>*h*</sub><sup>2</sup> at day *d* − 1, and *u*<sub>*</sub> and *B*<sub>0</sub> at day *d*. The reason is that *h* and
*N*<sub>*h*</sub><sup>2</sup> represent an initial state, with a MLD *h* in which thermals can develop underlying a stable stratification *N*<sub>*h*</sub><sup>2</sup>.
This initial state is modified by a whole day of forcing of *u*<sub>*</sub> and *B*<sub>0</sub>. Hence, the temporal evolutions of the MLD
∂<sub>t</sub>*h* are computed over this 1-day period (day *d* − 1 to day *d*). Finally, we note that the calculation of Δρ = ρ
(z = −1.1 hr) − ρ(z = −*h*), needed to obtain *N*<sub>*h*</sub><sup>2</sup> (see Equation 12), can present two problems. First, *N*<sub>*h*</sub><sup>2</sup> can be
negative if there is an instability at the base of the mixed layer. These cases represent less than 0.01 % of the points
and are simply discarded. Second, the calculation of Δρ is not defined if the mixed layer reaches the bottom of the
domain. These points, which correspond to a zero ∂<sub>t</sub>*h* evolution, also represent less than 0.01 % of all cases and are
discarded as well. As soon as these points are removed, the absolute sample size of the different simulations
(which depends on the horizontal resolution, the number of vertical levels and the total time of simulation) vary
from ≈130,000 points for the smallest one (the eORCA1 simulation extracted in the Western Mediterranean
region) to ≈20,000.000 for the largest one (the eORCA12 simulation extracted in the Labrador region).

## 2.9. Visualization in the Parameter Space

We characterize the MLD evolution through the relative change of the MLD ∂_t/h/h, expressed in %/day. The normalization by h allows us to compare different orders of magnitude of MLDs (which vary widely depending on different regions in the world, and/or different periods of the year). It would also have been possible to normalize with respect to time by dividing by a characteristic time t_c. Several possibilities were tested. However, since this normalization makes it more difficult to understand the variable, and since none of the trials produced any improvement in the results, none of the possibilities were retained.

The parameter space has three dimensions: λs, Rhy, and f/N. For exploring these three dimensions, we use pro-
jections into two-dimensional parameter spaces λs − f/N, Rh − f/N, and λs − Rh. To facilitate the intercomparison
of two graphs, we use hexagonal bin plots rather than scatter plots and we define four MLD evolution regimes
according to the value of ∂_t/h/h. Comparing two graphs can then be done by looking at the MLD evolution regime
obtained hexagon by hexagon. The four MLD evolution regimes are defined as follows

* ∂_t/h/h ≥ 10%/day: **Strong Deepening**

- 1%/day ≤ ∂_t/h/h < 10%/day: **Deepening**

* -1%/day ≤ ∂<sub>t</sub>/h/h < 1%/day: **Stable**

* $\partial_t h/h < -1\%/day$: **Restratification**

The regime of a hexagon is determined by the majority regime of its constituent points. That is, for every hexagon of a 2D parameter space we sort all the points inside this hexagon in one of the four regimes and the regime of the hexagon is the one that is the most represented. If this regime represents more than 75 % of the points, it is tagged as “highly significant.” If this percentage is between 50 % and 75 %, it is tagged as “significant.” If it is below 50 %, it is considered not significant. For statistical reasons, a hexagon is kept only if it contains at least 30 points.

**3. Results**

In this section, we populate the 3D parameter space (λs, R<sub>h</sub>, f/N<sub>h</sub>) with 1D simulations performed at the global scale and we show that considering only the two dimensions (λs, R<sub>h</sub>) is the best choice of 2D projection of this 3D parameter space. We then focus on this two-dimensional λs − R<sub>h</sub> parameter space and we present how this framework can be used to analyze 3D realistic ocean simulations. The first application is about the impact of the lateral resolution on the MLD evolution regimes. The second application focuses on the effect of the GM
parameterization which aims at representing the impact of the unresolved mesoscale processes in a coarse-resolution ocean model. Finally, a proof of concept of using observations as a truth is given in Appendix D.

### 3.1. Evaluation of the Three-Dimensional Parameter Space

The three-dimensional parameter space $\lambda_s - R_h - f/N_h$ is evaluated with 1D simulations performed at the locations of the 15° global grid (described in Section 2.7). We first show that the ($\lambda_s, R_h$) projection is the best choice of 2D projection of this 3D parameter space. We then highlight the influence of $f/N_h$.

Figure 2 displays the three two-dimensional projections of the 3D parameter space $\lambda_s - R_h - f/N_h$. Among these three projections, the 2D parameter space $\lambda_s - R_h$ exhibits the highest significance with the MLD evolution regimes of its hexagons being significant in 96 % of the cases. In comparison, the significance is 73 % for the $f/N_h - R_h$ projection and 73 % for the $\lambda_s - f/N_h$ projection (definition of the “significance” in Section 2.9). Moreover, among these significant hexagons, the $\lambda_s - R_h$ parameter space shows the highest number of highly significant hexagons: 74 % of all hexagons (23 % for $f/N_h - R_h$; 30 % for $\lambda_s - f/N_h$; cf. Table 2).

In addition to the high significance levels, the $\lambda_s - R_h$ projection also exhibits the best “spatial coherence”: the four MLD evolution regimes are organized in well-delimited continuous zones. The main thresholds delineating these zones are plotted in Figure 2a and are as follows:

* **Vertical threshold at $\lambda_s = -3$:** Physically, we expect restratification for $G/P|_{surf} < -1$, indicating that the surface buoyancy flux removes more TKE than the amount produced by the wind. In terms of $\lambda_s$ (i.e., $\propto G/P|_{surf}$), the threshold seems to be around $\lambda_s \approx -3$. The criterion $\lambda_s < -3$ corresponds to stable or restratification regimes and is consistent with a TKE-loss situation. This $\lambda_s < -3$ threshold is also observed in the $\lambda_s - f/N_h$ parameter space and so does not depend on $f/N_h$.

* **Horizontal thresholds in the range $-3 < \lambda_s < 0.2$:** In the range $-3 < \lambda_s < 0.2$, the boundaries are horizontal, indicating that, when the wind dominates over the buoyancy flux, only the value of $R_h$ is important for predicting the MLD evolution regime. A value of $R_h > 1,000$ corresponds to “stable” regime, $1,000 < R_h < 300$ corresponds to “deepening” regime and $R_h < 300$ corresponds to “strong deepening” regime. This progression according to $R_h$ corresponds to the traditional interpretation of a Richardson number. For high values of $R_h$, the shear $u_s/h$ is too weak to erode the pycnocline stratification $N_h$, leading to a stable regime. In contrast, low values of $R_h$ result in MLD deepening regimes.

* **Sloping lines in the zone $\lambda_s > 0.2$:** For $\lambda_s > 0.2$, demarcations follow lines with a 2/3 slope, which are isolines of $R_h^* = (N_h h/w^*)^2$. In this zone, only the value of the buoyancy-flux related Richardson number $R_h^*$ is hence important for predicting the MLD evolution regime. A value $R_h^* > 3000$ indicates a stable regime, while low values suggest deepening or strong deepening regimes. The fact that $R_h^*$ is the important dimensionless number in the $\lambda_s > 0.2$ zone informs us that this zone is a buoyancy-flux-dominant zone ($G > P$).

To summarize the previous points, $\lambda_s$ indicates a restratifying TKE-loss zone for $\lambda_s < -3$ and a TKE-gain zone for $\lambda_s > -3$. The TKE-loss zone is buoyancy-flux-dominant whereas the TKE-gain zone is either wind-dominant for $-3 < \lambda_s < 0.2$ and represented by demarcations by $R_h$, or buoyancy-flux-dominant for $\lambda_s > 0.2$ and represented by demarcations defined with $R_h^*$. For clarity, these interpretations based on $G$ and $P$ are added on the top of Figure 2a.

While the $\lambda_s - R_h$ projection is the best 2D projection of the 3D parameter space, we also explore the third dimension within this space. We sort the results according to their $f/N_h$ values and we plot different “slices” of the parameter space in Figure 3.

The influence of rotation, as assessed by the parameter $f/N_h$, appears to limit the deepening of the mixed layer. To illustrate this effect, we highlight in Figure 3 the demarcations between the stable regime and the deepening regime. We also plot in Figure 4 the value of these thresholds $R_{h,c}$ and $R_{h,c}^*$ as a function of $f/N_h$ (normalized by their values $R_{h,c0}$ and $R_{h,c0}^*$ for $f/N_h \in [10^{-3.5}; 10^{-3.0}]$). The stronger $f/N_h$, the lower are $R_{h,c}$ and $R_{h,c}^*$. Consequently, in the presence of rotation, a weaker stratification and/or a higher forcing ($u^*/h$ or $w^*/h$) are required to achieve the same level of deepening as without rotation. This can also be seen as a reduction of the effective impact of the forcing. For the wind, this reduced effective wind power input could be attributed to the generation of inertial oscillations that redistribute the horizontal momentum (Pollard et al., 1973). Concerning the buoyancy forcing, in such 1D simulations, rotation only appears in the mean horizontal velocities equations and therefore there is no physical mechanisms that could explained a reduced effective surface buoyancy flux (and indeed, free
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Results of the 1D simulations performed at the locations of the 15° global grid, plotted in the three different 2D projections of the 3D parameter space: (a) $\lambda_s - R_h$, (b) $f/N_h - R_h$, and (c) $\lambda_s - f/N_h$. The mixed layer depth (MLD) evolution regimes are defined based on the values of $\partial_t h/h$. The “strong deepening” regime is defined by $\partial_t h/h \ge 10\%/day$, “deepening” by $1\%/day \le \partial_t h/h < 10\%/day$, “stable” by $-1\%/day \le \partial_t h/h < 1\%/day$ and “restratification” by $\partial_t h/h < -1\%/day$. The regime of a hexagon is the majority regime of its constituent points. If this majority regime represents less than 50 % of the constituent points, the hexagon is tagged “non-significant” (superimposed black hexagon), if it is between 50 % and 75 % it is “significant” (superimposed black dot) and if it is higher than 75 % it is “highly significant” (nothing superimposed). A grid representing the slope 2/3 isolines of $R_h^*$ is added in the $\lambda_s > 0$ panel of the $\lambda_s - R_h$ parameter space. Dashed lines highlight demarcations between MLD evolution regimes discussed in the text.</figcaption></figure>
<table><caption>**Table 2**
Significances of the Mixed Layer Depth Evolution Regimes of the Hexagons for the Three 2D Projections of the 3D Parameter Space $\lambda_s - R_h - f/N_h$</caption><thead><tr><th>2D parameter space</th><th>λ<sub>s</sub> - R<sub>h</sub> (%)</th><th>f/N<sub>h</sub> - R<sub>h</sub> (%)</th><th>λ<sub>s</sub> - f/N<sub>h</sub> (%)</th></tr></thead><tbody><tr><td>Highly significant hexagons</td><td>74</td><td>23</td><td>30</td></tr><tr><td>Significant hexagons</td><td>22</td><td>50</td><td>43</td></tr><tr><td>Not significant hexagons</td><td>4</td><td>27</td><td>27</td></tr></tbody></table>

convective simulations performed with or without rotation give exactly the same results). Hence, the reduction of the deepening in the buoyancy-flux-dominant zone must be also attributed to a reduction of the wind power input even if it is not the dominant forcing.

As the effect of rotation can be seen as a reduction of the effective wind power input, its impact is more pronounced in the wind-dominant zone and hence $R_{h,c}$ decreases more with $f/N_h$ than $R_{h,c}^*$ (cf. Figure 4). Therefore, the region for which the MLD deepening is driven by the wind narrows with $f/N_h$ compared to the region for which the MLD deepening is driven by the surface buoyancy

forcing. The $\lambda_{s,c}$ thresholds that delineate these two regions are plotted in Figure 3 and their corresponding values are reported in Figure 4 (following Equation 11, it could be calculated by $\lambda_{s,c} = (R_{h,c}/R_{h,c}^*)^{3/2}$). This dependency on $f/N_h$ suggests that regions near the Equator are more likely to be in the wind-dominant regime, whereas high-latitude regions are more inclined toward a buoyancy-dominant regime.

As previously observed, the $\lambda_s - R_h$ parameter space (Figure 2a) exhibits high levels of significance even if the dimension $f/N_h$ is not considered. This suggests that variations of the $f/N_h$ parameter are less important than variations of $\lambda_s$ and $R_h$ for predicting a MLD evolution regime. Compared to the influence of $f/N_h$ observed in Figure 3, this behavior can be explained by the fact that the $f/N_h$ distribution of the 1D simulations at the global scale is not uniform but dominated by values of $f/N_h \in [10^{-2.5}; 10^{-1.5}]$ (cf. Figure 5; and one can indeed see that demarcations of Figure 2a are close to the demarcations in Figures 3c and 3d). Given the high statistical performance of the $\lambda_s - R_h$ parameter space, and because it is easier to work in two dimensions, we will focus solely on the $\lambda_s - R_h$ projection in the remainder of this article. The dimension $f/N_h$ will only be considered if necessary to comprehend low significance levels.

We conclude this section by noting that the statistical performance of the parameter space is not specific to the TKE vertical mixing scheme (see Appendix A for a brief presentation of results with the KPP scheme). Additionally, for informative purposes, the density maps and associated joint Probability Density Functions (PDF) showing the density distribution of the values of $\lambda_s$, $R_h$, and $f/N_h$ in the three 2D projections of the 3D parameter space are given in Appendix B. This information can be useful when selecting relevant values of forcing and preconditioning conditions ($u_*$, $B_0$, $N_h$) in the context of parameter tuning (Souza et al., 2020; Wagner et al., 2023). Finally, it should be noted that all the results of this section were obtained with 1D simulations. The equivalent diagnostics were also realized with a 3D model (the eORCA12 simulation described in Section 2.5) and are given in Appendix C. To summarize quickly here, the results about the influence of the $f/N_h$ parameter are hardly observable with eORCA12 because of a higher number of non-significant hexagons than the 1D simulations. However, it is confirmed that the $\lambda_s - R_h$ parameter space is the best choice of 2D projections of the 3D parameter space also for eORCA12 (cf. significance levels in Table C1).

## 3.2. Influence of the Horizontal Resolution on the MLD Evolution Regimes for 3D Ocean Circulation Models

Figure 6 displays the results of the 1D simulations, eORCA1, eORCA12, and eNATL60 in the $\lambda_s - R_h$ parameter space for the global scale and in the Western Mediterranean region. Since all four simulations used the same 1D vertical scheme TKE + EVD, any variations between the figures are attributed to the influence of lateral processes. We will not comment on the extent of the parameter space that is explored by the different simulations but rather comment on the differences of MLD evolution regimes observed in the regions explored in common. The reason is that the extent of the explored region of the parameter is dependent on the absolute sample size. For instance, the eORCA12 and eNATL60 simulations explore a larger region of the parameter space than the eORCA1 and 1D simulations in the Western Mediterranean case just because the sample size is bigger (longer time period and/or higher vertical resolution and/or higher number of vertical levels).

At the global scale, the main demarcation lines are consistent across the three simulations. This suggests that, for predicting the MLD evolution, lateral processes are of second importance in comparison with the 1D processes presented via $\lambda_s$ and $R_h$. However, this observation may not be locally valid. Extractions of the three same simulations in the Western Mediterranean, in addition to the eNATL60 simulation (which could not be considered at the global scale due to its basin-scale nature), reveal significant variations across different resolutions. This
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results of the 1D simulations performed at the locations of the 15° global grid plotted in the $\lambda_s - R_h$ parameter space. The results are filtered according to their $f/N_h$ values: (a) $f/N_h \in [10^{-3.5}; 10^{-3.0}]$, (b) $f/N_h \in [10^{-3.0}; 10^{-2.5}]$, (c) $f/N_h \in [10^{-2.5}; 10^{-2.0}]$, (d) $f/N_h \in [10^{-2.0}; 10^{-1.5}]$, and (e) $f/N_h \in [10^{-1.5}; 10^{-1.0}]$. Dashed lines highlight the thresholds $R_{h,c}$ and $R_{h,c}^*$ between the stable and the deepening regimes, and $\lambda_{s,c}$ the limit between the wind-dominant and the surface-buoyancy-flux-dominant zones. Values $R_{h,c0}$, $R_{h,c0}^*$ and $\lambda_{s,c0}$ are the ones for $f/N_h \in [10^{-3.5}; 10^{-3.0}]$. Other graphical conventions as in Figure 2.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dependence on $f/N_h$ of the three demarcation thresholds $R_{h,c}$, $R_{h,c}^*$ and $\lambda_{s,c}$. The thresholds $R_{h,c}$ and $R_{h,c}^*$ indicate the demarcation between the stable and the deepening regimes respectively in the wind-dominant zone and in the surface-buoyancy-forcing-dominant zone. The threshold $\lambda_{s,c} = (R_{h,c}/R_{h,c}^*)^{3/2}$ indicates the transition between these wind-dominant zone and surface-buoyancy-forcing-dominant zone. These three thresholds are plotted normalized by their values at the lowest $f/N_h$: $R_{h,c0}$, $R_{h,c0}^*$ and $\lambda_{s,c0}$. Letters (a), (b), (c), (d), and (e) refer to the subfigures of Figure 3.</figcaption></figure>

indicates that the lateral processes play a substantial role in the Western Mediterranean sea and cannot be neglected when compared to the 1D processes. The next four paragraphs detail these “significant variations” that are observed between the four simulations in the Western Mediterranean region.

For $\lambda_s > -3$ and high values $R_h > 4,000$ and $R_h^* > 2000$, the high resolution simulations (1/12° and 1/60°) mainly result in a restratification regime. High values of $R_h$ and $R_h^*$ indicate a stable MLD in terms of the 1D processes. The almost complete replacement of the stable regime by the restratification regime suggests that the lateral processes of restratification, such as restratification by baroclinic instability at convective fronts or at mesoscale eddy fronts, are predominant in this region and easily develop at timescales $\le 1$ day. In the same conditions ($\lambda_s > -3$; $R_h > 4,000$; $R_h^* > 2000$), the 1° simulation behaves similarly to the 1D simulations, exhibiting a “stable” regime. This suggests that the coarse-resolution 1° simulation poorly resolves these lateral processes of restratification.

For $\lambda_s > -3$ and low $R_h < 400$, going from 1D to 3D does not have a significant impact: the “strong deepening” regime is maintained in all four simulations. Low values of $R_h$ indicate an unstable water column where the stratification is low compared to the wind forcing. The preservation of the “strong deepening” regime in all simulations suggests that for $R_h < 400$ the lateral processes of restratification cannot neutralize this instability and, therefore, play a secondary role.

The zone with $\lambda_s < -3$, dominated by surface buoyancy fluxes is a zone of restratification or stable regimes for the 1D simulations. As mentioned earlier, we associated this behavior with a TKE-loss in the 1D TKE budget. Interestingly, this $\lambda_s < -3$ zone characterized by restratification or stable regimes is still observed in the 3D models. This implies that for a dominant surface-buoyancy restratifying flux ($\lambda_s < -3$), the lateral processes of TKE generation (such as an Ekman flow creating an equivalent destabilizing wind-driven buoyancy flux, see for example Thomas & Lee, 2005) are of secondary importance compared to the processes of the 1D TKE budget.

The percentage of significant hexagons decreases when the resolution increases: it is 97 % in 1D, 88 % at 1°, 84 % at 1/12° and 78 % at 1/60°. Non-significant hexagons indicate that the performance of the parameter space for

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Relative frequency of the values of $f/N_h$ for the 1D simulations at the global scale. The five slices used in Figure 3 are highlighted in light blue and the percentages of values falling in each of them are given.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Results in the $\lambda_s - R_h$ parameter space of the (a) 1D, (b) 1° eORCA1 and (c) 1/12° eORCA12 simulations at the global scale; and of the (d) 1D, (e) 1° eORCA1, (f) 1/12° eORCA12 and (g) 1/60° eNATL60 simulations for the Western Mediterranean region. Dashed lines highlight demarcations discussed in the text. Graphical conventions as in Figure 2.</figcaption></figure>

assessing the MLD evolution is hampered by the importance of lateral processes. Considering the parameter $f/N_h$ does not improve the results (not shown). Other parameters, some of which are described in Section 2.2, could improve the predictability, but we have not studied these higher-dimensional parameter spaces further here.

### 3.3. Impact of the GM Parameterization on a 1° Coarse-Resolution Simulation

In Figure 7, we plot the results of the eORCA1, eORCA1GM, and eNATL60 simulations in the Western Mediterranean and the Labrador regions separately. We recall that eORCA1GM differs from eORCA1 solely due to the addition of the GM parameterization, designed to represent the adiabatic advective effect of unresolved mesoscale processes (Gent, 2011). Using eNATL60 as a reference helps evaluate how these mesoscale processes can influence the MLD evolution regimes in the parameter space. Therefore, comparing eORCA1GM and eNATL60 provides valuable insights into the impact of the GM parameterization, even though the GM parameterization was not designed to tackle the impact of the mesoscale processes in the mixed layer.

The impact of the GM parameterization is minimal in the Western Mediterranean region (Figures 7a and 7b), although resolving the mesoscale processes changes for instance the MLD evolution regimes from “stable” to “restratification” at high $R_h > 2,000$ and high $R_h^* > 2000$ (Figures 7a and 7c).
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Results in the Western Mediterranean region of the (a) eORCA1, (b) eORCA1GM, and (c) eNATL60 simulations in the $\lambda_s - R_h$ parameter space. Panels (d–f) are the same in the Labrador region. Dashed lines highlight important demarcations discussed in the text. Graphical conventions as in Figure 2.</figcaption></figure>

In the Labrador region, the GM parameterization has a visible impact where its main contribution is to generate a restratification zone for intermediate stability conditions ($200 < R_h < 2,000$ in $-3 < \lambda_s < 2$; $200 < R_h^* < 2000$ in $\lambda_s > 2$) whereas it was before mainly a deepening zone (Figures 7d and 7e). However, neither the stable zone ($R_h > 2,000$ in $-3 < \lambda_s < 2$; $R_h^* > 2000$ in $\lambda_s > 2$) nor the strong deepening zone ($R_h < 200$ in $-3 < \lambda_s < 2$; $R_h^* < 200$ in $\lambda_s > 2$) are affected. These changes are not comparable with the effect of the mesoscale processes represented by the 1/60° results (comparison Figures 7d and 7f) for which, for instance, the restratification zone at intermediate stability conditions ($200 < R_h < 2,000$ in $-3 < \lambda_s < 2$; $200 < R_h^* < 2000$ in $\lambda_s > 2$) is not observed.

Therefore, these two cases highlight that the impact of mesoscale processes on the MLD evolution is probably not adequately captured by GM. To better characterize the impact of the GM parameterization as a function of the position in the parameter space, we can examine a proxy for its activation. The GM parameterization tends to flatten isopycnals by advecting tracers via eddy-induced velocities (Gent et al., 1995)
$$ 
\begin{aligned}
u_{GM} &= -\partial_z(\kappa_{GM}S_x) \\
v_{GM} &= -\partial_z(\kappa_{GM}S_y) \\
w_{GM} &= \partial_x(\kappa_{GM}S_x) + \partial_y(\kappa_{GM}S_y)
\end{aligned} 
\qquad (13) $$

with $\kappa_{GM}$ the isopycnal thickness diffusivity, $S_x = -\partial_x\rho/\partial_z\rho$ the zonal isopycnal slope and $S_y = -\partial_y\rho/\partial_z\rho$ the meridional isopycnal slope.

We construct an index to quantify the magnitude of the GM rectification by considering the horizontal transports integrated over the mixed layer $\gamma_x = \int_{-h}^{0} -\partial_z(\kappa_{GM}S_x) dz$ and $\gamma_y = \int_{-h}^{0} -\partial_z(\kappa_{GM}S_y) dz$. The surface boundary condition imposes $w_{GM} = 0$. This condition is often satisfied by taking $\kappa_{GM}S_x = \kappa_{GM}S_y = 0$ at the surface. Thus $\gamma_x = \kappa_{GM}(z = -h)S_x(z = -h)$ and $\gamma_y = \kappa_{GM}(z = -h)S_y(z = -h)$: the integrated horizontal transports are proportional to the isopycnal slopes at the mixed layer base. The index is then constructed as the maximal isopycnal slope over the x and the y axes

$$ 
\begin{aligned}
S_h &= \max(|S_x(z = -h)|; |S_y(z = -h)|) \\
&= \max\left(\left|\frac{\partial_x \rho(z = -h)}{\partial_z \rho(z = -h)}\right|; \left|\frac{\partial_y \rho(z = -h)}{\partial_z \rho(z = -h)}\right|\right)
\end{aligned} 
\qquad (14) $$

with $\partial_z\rho(z = -h)$ calculated over a distance 0.1 hr below the mixed layer: $\partial_z\rho(z = -h) = \rho(z=-1.1 \frac{h}{0.1h}) - \rho(z=-h)$.

In the Labrador region, the previously-described zone ($200 < R_h < 2,000$ in $-3 < \lambda_s < 2$; $200 < R_h^* < 2000$ in $\lambda_s > 2$) where the GM parameterization has a notable impact is characterized by values of $S_h > 0.5$ m/km (Figures 8e, 8f, and 8h). We therefore decided to define a threshold $S_h = 0.5$ m/km below which the GM parameterization is expected to have a negligible impact. This is supported by the fact that the values of $S_h$ found in the Western Mediterranean region are mainly below this threshold (Figure 8a). Filtering the eORCA1GM results by the condition $S_h < 0.5$ m/km yields results comparable to the results of the eORCA1 simulation (comparison Figures 8c and 8d, and Figures 8g and 8h). This confirms that the observed impacts of the GM parameterization can be understood through the values of $S_h$. Hence, this number is an additional dimensionless number that could be considered when $R_h$ and $\lambda_s$ alone do not provide robust predictions. In future studies, exploring projections into parameter spaces ($\lambda_s, S_h$) or ($R_h, S_h$) could be promising avenues.

## 4. Conclusions and Discussion

This study introduces a three-dimensional parameter space designed to facilitate the analysis and the intercomparison of the ocean MLD evolution between numerical models. The parameter space consists of three dimensionless numbers $R_h$, $\lambda_s$, and $f/N_h$ derived through dimensional analysis: $\lambda_s$ evaluates the relative influence of the buoyancy forcing and the wind forcing for producing/destroying TKE in the surface layer, $R_h$ is the Richardson number describing the competition between the stabilizing effect of the pycnocline stratification and the destabilizing impact of the wind-induced shear. Finally, $f/N_h$ evaluates the influence of the Earth's rotation.

The $\lambda_s - R_h - f/N_h$ parameter space was first evaluated in the context of 1D simulations. Four MLD evolution regimes were defined based on the value of the relative MLD change $\partial_h/h$: “Restratification,” “Stable,” “Deepening,” and “Strong Deepening.” We showed that the influence of rotation tends to stabilize the water column by reducing the effective forcing of the wind. Consequently, MLD deepening in high latitude regions is more inclined to be dominated by the surface buoyancy forcing whereas the Equator is more inclined to be wind-dominant. This can be related to the numerous studies that assess the relative importance of wind and surface buoyancy forcings in different regions (Dong et al., 2007; Downes et al., 2011; Gao et al., 2023; Holte et al., 2012; Sallée et al., 2010, 2021). For example, Sallée et al. (2010) stated that the surface buoyancy forcing in the Southern Ocean (high latitude) dominates the wind forcing by one order of magnitude. Our study, which shows the dependence on $f/N_h$ of the $\lambda_{s,c}$ threshold between the two regimes, provides a new practical way of determining the relative importance of surface buoyancy flux versus wind.

The influence of the $f/N_h$ parameter is less important than the $\lambda_s$ and $R_h$ parameters for predicting the MLD evolution. The 2D parameter space $\lambda_s - R_h$ indeed exhibits high statistical performances with, in 96 % of the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Results in the Western Mediterranean region in the $\lambda_s - R_h$ parameter space of (a) $S_h$ in eORCA1, (b) mixed layer depth (MLD) evolution regimes in eORCA1GM, (c) MLD evolution regimes in eORCA1GM restricted to $S_h < 0.5$ m/km and (d) MLD evolution regimes in eORCA1. Panels (e–h) are the same in the Labrador region. The slope $S_h$ of the isopycnes at the MLD is expressed in m/km. In panel (e), values that exceeded 1.5 m/km are represented in green. Dashed lines highlight demarcations discussed in the text. Graphical conventions as in Figure 2.</figcaption></figure>

cases, a pair $(\lambda_s, R_h)$ that corresponds to a unique MLD evolution regime. In other words, instead of examining all the preconditioning and forcing conditions, one can just calculate the two dimensionless numbers $\lambda_s$ and $R_h$ for predicting the MLD evolution regimes. Also, the MLD evolution regimes appear in well-separated zones. This spatial coherence of regimes in this parameter space allows us to define thresholds on $(\lambda_s, R_h)$ to predict MLD evolution regimes.

The thresholds of the $\lambda_s - R_h$ parameter space were described in the context of 1D simulations. The criterion $\lambda_s < -3$ indicates stable or restratification regimes and is valid for all $f/N_h$ values. The wind-dominant zone $-3 < \lambda_s < 0.2$ is characterized by transitions according to $R_h$-only thresholds. In the buoyancy-flux-dominant zone $\lambda_s > 0.2$, transitions between regimes can be seen as thresholds on $R_h^* = (N_h h / w_*)^2$, the Richardson number associated with the destabilizing buoyancy flux. This threshold at $\lambda_s \approx 0.2$ between the wind-dominant and the surface-buoyancy-flux-dominant zones is the one for the global scale, which is representative of $f/N_h$ values in $[10^{-2.5}; 10^{-1.5}]$, and must be adjusted for different values of $f/N_h$.
Two applications of the parameter space were presented and we show how it may be used with realistic ocean models. In the first application, we intercompare ocean simulations at different horizontal resolutions to evaluate the effect of lateral processes on the MLD evolution. We showed that lateral processes play a secondary role for low values of $R_h$ and $R_h^*$: the stratification effect is weak compared to the forcing $u_*$ or $w_*$ and the “strong deepening” regime can be predicted without considering lateral processes. When the water column is stable with respect to 1D processes (large values of $R_h$ and $R_h^*$), we showed that the lateral restratification processes become dominant in the high-resolution simulations (1/12° and 1/60°). These lateral processes may for example, include the restratification by baroclinic instability at convective fronts or at mesoscale eddy fronts. However, the 1° simulation behaves as the 1D model, suggesting that the lateral processes of restratification are not resolved at this resolution. Finally, in the high-resolution simulations (1/12° and 1/60°), the non-significant zones at mid values of $R_h$ and $R_h^*$ indicate that the lateral processes are dominant and that other dimensionless numbers could be considered for predicting the MLD evolution regime.

The second application shows that the adiabatic advective effect of the mesoscale processes parameterized by GM parameterization does not capture the full impact of unresolved mesoscale processes on the MLD evolution regimes in a coarse-resolution 1° simulation. In this context, we introduced the dimensionless number $S_h$ which is the maximal isopycnal slope at the mixed layer base. This slope is one of the other dimensionless numbers that could be considered when the two ($R_h$, $\lambda_s$) are not sufficient for obtaining robust predictions. Particularly, projections into the parameter space ($\lambda_s$, $S_h$) or ($R_h$, $S_h$) could constitute some developments for future works.

The two applications presented in this study are not exhaustive. We decided to focus here on the use of the parameter space for model sensitivity studies. Future work could use the parameter space for comparing the behaviors of different vertical mixing schemes (KPP, TKE, GLS) and for comparing coupled and forced models. The information of the joint PDF of the three 2D projections of the 3D parameter space, given in Appendix B could also be used for choosing relevant values of forcing and preconditioning conditions ($u_*$, $B_0$, $N_h$) in the context of parameter tuning (Souza et al., 2020; Wagner et al., 2023). Beyond these direct applications, an interesting extension of the approach would be to evaluate the performance of the parameter space with LES data and observations. In Appendix D, we propose a short proof of concept using observations from the PAPA station. In the future, we could greatly increase the sample size of our observation data by using, for instance, ARGO float data. These floats give profiles over 10-day periods and the parameter space will need to be assessed with this new period. For the LES, it would be possible to keep the 1-day period developed in this study or to try also with shorter or longer periods. If the statistical performance is still obtained in these contexts, the parameter space could become an informative tool for calibrating the mixing schemes using LES or observational data as a truth.

Last, we want to mention the study of Johnson et al. (2023) that proposes an interesting methodology to evaluate the robustness and uncertainty of 1D models. They describe a framework to highlight situations where different models will predict a wide spread in the evolution of the mixed layer. While they applied their methodology to a limited event (20 days at a specific location), we envision using the same methodology but for all the cases described in the present article. The metrics they define can be visualized in our parameter space, and we would then see the emergence of “uncertain” zones where the definition of “uncertainty,” thanks to their method, could be directly linked to either model formulation, model state, or model parameters. This research avenue is currently being investigated in our research group.

## Appendix A: Analysis With the KPP Vertical Mixing Scheme

To verify if the statistical performance of the parameter space is sensitive to the vertical mixing scheme, we performed the same collection of 1D water column simulations described in Section 2.7 but with a KPP scheme instead of a TKE + EVD scheme. Figure A1 presents the results of these simulations with the same conventions as Figure 2. Again, the $\lambda_s - N_h$ parameter space performs well with 96 % of significant hexagons, and spatial coherence of well-delimited zones is still obtained. The demarcation thresholds (represented by dashed lines) could again be discussed. In short, all the diagnostics we have done previously could have been done with simulations based on the KPP scheme as well, and future research could focus on analyzing the difference in behaviors between the TKE + EVD scheme and the KPP mixing scheme. For example, by comparing Figures 2 and A1, we can see that the demarcations between the “strong deepening” and the “deepening” regimes in the $-3 < \lambda_s < 0.2$ region is, with
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Same as Figure 2 but with 1D simulations using a K-Profile Parameterization vertical mixing scheme instead of a Turbulent Kinetic Energy + Enhanced Vertical Diffusivity scheme.</figcaption></figure>

the KPP scheme, close to $R_h \approx 400-500$ rather than $R_h \approx 300$ with the TKE + EVD scheme. In the wind-dominant region, at mid-$R_h$ values, the KPP scheme is hence more inclined to deepen the MLD than the TKE + EVD scheme.

## Appendix B: Joint PDF of Three 2D Projections of the 3D Parameter Space

We plot in Figure B1 the density maps in the three 2D projections of the 3D parameter space and the contours of the associated joint PDF, calculated with the 1D simulations outputs at the global scale through the Python functions provided by Q. Li et al. (2019). These contours enclose 30 % (black), 60 % (blue), 90 % (green), and 99 % (yellow) of all instances centered at the highest PDF.
<figure><img src="image_13.png" /><figcaption>**Figure B1.** Density maps of the 1D simulations at the global scale for the (a) $\lambda_s - R_h$, (b) $f/N_h - R_h$ and (c) $\lambda_s - f/N_h$ parameter space. The contours of the associated joint Probability Density Functions (PDF) are superimposed. These contours enclose 30 % (black), 60 % (blue), 90 % (green), and 99 % (yellow) of all instances centered at the highest PDF.</figcaption></figure>
## Appendix C: Evaluation of the Three-Dimensional Parameter Space With eORCA12

All the results presented in Section 3.1 were realized with the 1D simulations. It is interesting to know if we can obtain the same results with a 3D model (we use the eORCA12 simulation described in Section 2.5). Figure C1, Table C1, and Figure C2 present the same diagnostics as Figure 2, Table 2, and Figure 2, but with eORCA12 rather than the 1D simulations. For indicative purpose, the demarcation thresholds found for the 1D case are again plotted in the figures.

<figure><img src="image_14.png" /><figcaption>**Figure C1.** Same as Figure 2 but for the eORCA12 simulation rather than the 1D simulations. For indicative purpose, the green dashed lines highlight the same demarcations as the ones represented in Figure 2.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure C2.** Same as Figure 3 but for the eORCA12 simulation rather than the 1D simulations. For indicative purpose, the green dashed lines highlight the same demarcations as the ones represented in Figure 3.</figcaption></figure>
<caption>Table C1</caption>

Significances of the Mixed Layer Depth Evolution Regimes of the Hexagons for the Three 2D Projections of the 3D Parameter Space $\lambda_s - R_h - f/N_h$ (Same as Table 2 but for the eORCA12 Simulation Rather Than the 1D Simulations)

<table><thead><tr><th>2D parameter space</th><th>λ<sub>s</sub> - R<sub>h</sub> (%)</th><th>f/N<sub>h</sub> - R<sub>h</sub> (%)</th><th>λ<sub>s</sub> - f/N<sub>h</sub> (%)</th></tr></thead><tbody><tr><td>Highly significant hexagons</td><td>41</td><td>16</td><td>20</td></tr><tr><td>Significant hexagons</td><td>42</td><td>37</td><td>27</td></tr><tr><td>Not significant hexagons</td><td>17</td><td>47</td><td>53</td></tr></tbody></table>

With Figure C1 and Table C1, we again conclude that the $\lambda_s - R_h$ parameter space is the best choice of 2D projections of the 3D parameter space. This projection indeed presents the lowest percentage of non-significant hexagons and the highest spatial coherence with well-delimited zones of distinct MLD evolution regimes.

On the other hand, Figure C2 does not allow us to generalize the results found with the 1D simulations. Indeed, because of the higher number of non-significant hexagons, the clear demarcations between the “deepening” and the “stable” regimes can be identified only in Figures C2c and C2d, and the dependence on $f/N_h$ of the three demarcation thresholds $R_{h,c}$, $R_{h,c}^*$ and $\lambda_{s,c}$ can hence not be drawn for eORCA12 (as we did in Figure 4 for the 1D simulations). Realizing new extractions of a 3D simulation with more data points (either with a finer spatial resolution or with a larger time coverage) would allow us to conclude on these points. For this study, we did not need it as we decided to focus the $\lambda_s - R_h$ parameter space to assess and compare 3D simulations (Sections 3.2 and 3.3).

## Appendix D: Mixed Layer Dynamics at the PAPA Station

In this appendix, we highlight the possibility of using the parameter space with observations. We focus here on a subset of the data gathered at the PAPA station. It consists of 11 years (from 2010 to 2020) of temperature and salinity vertical profiles along with the surface heat flux, fresh water flux and wind stress. This test case is one of those available in the General Ocean Turbulence Model (GOTM, Burchard et al., 1999; Umlauf & Burchard, 2005, www.gotm.net). More information about how bulk formulae for the surface fluxes have been used is discussed in Burchard et al. (1999). The goal here is to use these observations as a truth for evaluating simulation results. Therefore, we selected the closest grid point to the PAPA station in the eORCA12 simulation and we extracted the associated data for the full 10-year period (2006–2015).

In Figure D1, we plot the MLD evolution regimes in the $\lambda_s - R_h$ parameter space both for the PAPA observations and for the eORCA12 simulation. Even if the number of non-significant hexagons is higher for the observations than the model (31 % vs. 6 %), the significant hexagons for the observations still exhibit a “spatial coherence” with well delimited zones. Considered as a truth, these demarcations that we discuss in the next paragraphs give insight on how the model performs for representing the MLD evolution.

For $\lambda_s < -3$, the observations result in a restratification regime. This behavior of giving non-deepening regimes for $\lambda_s < -3$ is also observed for the model that results in the stable or the restratification regimes. Taking a step back, this behavior is even shared by all the model results shown in this paper (with an exception for the eORCA1 model that presents a small deepening band at low $R_h$). This threshold $\lambda_s = -3$ is therefore robust and could be considered in the future as a predicting tool.

In unstable conditions ($\lambda_s > -3$ and $R_h < 600$), the model performs well compared to the observations, with both exhibiting deepening regimes. The only difference is about the intensity of this deepening which is a bit higher in the observations where only the strong deepening regime is obtained whereas the model gives a deepening regime close to the $R_h = 600$ demarcation.

For high-stability conditions ($\lambda_s > -3$ and $R_h > 7,000$), the model also performs well by giving the stable regime such as the observations.

The major differences between the model and the observations are located at intermediate stability conditions ($\lambda_s > -3$ and $600 < R_h < 7,000$). In these conditions, the model mainly results in stable and deepening regimes while the observations mostly indicate a restratification regime. We postulate that this restratification is mainly due to lateral processes (at least in the $\lambda_s > 0$ panel where no 1D process can explain a restratification) and the
<figure><img src="image_16.png" /><figcaption>**Figure D1.** Comparison in the $\lambda_s - R_h$ parameter space of (a) observations at the PAPA station and (b) eORCA12 simulation extracted at the closest grid point to this location. Dashed lines highlight demarcations discussed in the text. Graphical conventions as in Figure 2.</figcaption></figure>

intensity of these lateral processes are hence not well captured in the model. More generally, it is interesting to note that the restratification occurs in these intermediate stability conditions rather than in more stable water columns. This non-trivial behavior was already observed in some of the simulations results (see Figures 6f, 6g, 7e, and 7f).

This first exploration provides a proof of concept for the use in the parameter space of observations as a truth. In the future, we will work on greatly increasing the number of observations and LES data that we will use. We hope to provide interesting insights on how the different models perform for representing the MLD evolution.

## Data Availability Statement

All the codes used for the study are published on GitHub (Legay, 2024b, https://github.com/legaya/James2024ParameterSpace). This includes the Jupyter Notebook used for performing the 1D simulations and all the analyses, the 1D model described in Section 2.1 as a Fortran Module "scm_oce.so," and the Fortran codes needed for generating this module. The observations at the PAPA station, as well as the eNATL60, eORCA12, eORCA1GM, eORCA1, and 1D simulations outputs needed for realizing the figures are available as netCDF files and "npz" archives (Legay, 2024a).
Acknowledgments

We sincerely thank Brandon Reichl, Mike
Bell and three anonymous reviewers for
their valuable time and helpful feedback.
Their input has been instrumental in
enhancing this manuscript. We thank
Florian Lemarié for providing the 1D code,
and Aurélie Albert for the extractions of
the eNATL60 and eORCA12 simulations.
All the computations presented in this
paper were performed using the GRICAD
infrastructure (https://gricad.univ-
grenoble-alpes.fr), which is supported by
Grenoble research communities. Thierry
Penduff is supported by the MixED Layer
hEterogeneitY (MEDLEY) project which
has received funding from the Joint
Programming Initiative (JPI) Climate and
JPI Oceans programs under the 2019 joint
call, managed by the French Agence
Nationale de la Recherche (contract no. 19-
JPOC-0001-01).

References

Aurnou, J. M., Horn, S., & Julien, K. (2020). Connections between nonrotating, slowly rotating, and rapidly rotating turbulent convection transport scalings. *Physical Review Research, 2*(4), 043115. https://doi.org/10.1103/PhysRevResearch.2.043115

Banks, H. T., & Gregory, J. M. (2006). Mechanisms of ocean heat uptake in a coupled climate model and the implications for tracer based predictions of ocean heat uptake. *Geophysical Research Letters, 33*(7), L7608. https://doi.org/10.1029/2005GL025352

Belcher, S. E., Grant, A. L. M., Hanley, K. E., Fox-Kemper, B., Van Roekel, L., Sullivan, P. P., et al. (2012). A global perspective on Langmuir turbulence in the ocean surface boundary layer: FRONTIER. *Geophysical Research Letters, 39*(18), L18605. https://doi.org/10.1029/2012GL052932

Bernardello, R., Marinov, I., Palter, J. B., Galbraith, E. D., & Sarmiento, J. L. (2014). Impact of Weddell Sea deep convection on natural and anthropogenic carbon in a climate model. *Geophysical Research Letters, 41*(20), 7262–7269. https://doi.org/10.1002/2014GL061313

Boccaletti, G., Ferrari, R., & Fox-Kemper, B. (2007). Mixed Layer instabilities and restratification. *Journal of Physical Oceanography, 37*(9), 2228–2250. https://doi.org/10.1175/JPO3101.1

Bouillaut, V., Lepot, S., Aumartre, S., & Gallet, B. (2019). Transition to the ultimate regime in a radiatively driven convection experiment. *Journal of Fluid Mechanics, 861*, R5. https://doi.org/10.1017/jfm.2018.972

Brodeau, L., Le Sommer, J., & Albert, A. (2020). Ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. *Zenodo*. https://doi.org/10.5281/ZENODO.4032732

Burchard, H., & Bolding, K. (2001). Comparative analysis of four second-moment turbulence closure models for the oceanic mixed layer. *Journal of Physical Oceanography, 31*(8), 1943–1968. https://doi.org/10.1175/1520-0485(2001)031(1943:CAOFSM)2.0:CO;2

Burchard, H., Bolding, K., Villarreal, M., Centre, E. C. J. R., & Institute, S. A. (1999). GOTM, a general ocean turbulence model: Theory, implementation and test cases. Space Applications Institute. Retrieved from https://books.google.fr/books?id=zsJUHAAACAAJ

Damerell, G. M., Heywood, K. J., Calvert, D., Grant, A. L., Bell, M. J., & Belcher, S. E. (2020). A comparison of five surface mixed layer models with a year of observations in the North Atlantic. *Progress in Oceanography, 187*, 102316. https://doi.org/10.1016/j.pcean.2020.102316

de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A., & Iudicone, D. (2004). Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. *Journal of Geophysical Research, 109*(C12), C12003. https://doi.org/10.1029/2004JC002378

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., et al. (2011). The ERA-interim reanalysis: Configuration and performance of the data assimilation system. *Quarterly Journal of the Royal Meteorological Society, 137*(656), 553–597. https://doi.org/10.1002/qj.828

Deremble, B. (2016). Convective plumes in rotating systems. *Journal of Fluid Mechanics, 799*, 27–55. https://doi.org/10.1017/jfm.2016.348

Dingwall, J., & Taylor, J. R. (2024). Model for the cyclonic bias of convective vortices in a rotating system. *Physical Review Fluids, 9*(3), 033503. https://doi.org/10.1103/PhysRevFluids.9.033503

Dong, S., Gille, S. T., & Sprintall, J. (2007). An assessment of the Southern Ocean mixed layer heat budget. *Journal of Climate, 20*(17), 4425–4442. https://doi.org/10.1175/JCLI4259.1

Downes, S. M., Budnick, A. S., Sarmiento, J. L., & Farneti, R. (2011). Impacts of wind stress on the Antarctic circumpolar current fronts and associated subduction: Acc fronts and subduction. *Geophysical Research Letters, 38*(11), L11605. https://doi.org/10.1029/2011GL047668

Ekman, V. W. (1905). On the influence of Earth's rotation on ocean currents. *Arxiv for Matematik, Astronomi Och Fysik, 2*, 1–53.

Fearon, G., Herbette, S., Veitch, J., Cambon, G., Lucas, A. J., Lemarié, F., & Vichi, M. (2020). Enhanced vertical mixing in coastal upwelling systems driven by diurnal-inertial resonance: Numerical experiments. *Journal of Geophysical Research: Oceans, 125*(9), e2020JC016208. https://doi.org/10.1029/2020JC016208

Fofonoff, N. P., & Millard, R., Jr. (1983). *Algorithms for the computation of fundamental properties of seawater (technical report)*. UNESCO. https://doi.org/10.25607/OBP-1450

Fox-Kemper, B., Ferrari, R., & Hallberg, R. (2007). Parameterization of mixed layer eddies. I: Theory and diagnosis. *Journal of Physical Oceanography, 38*(6), 1145–1165. https://doi.org/10.1175/2007JPO3792.1

Frank, D., Landel, J. R., Dalziel, S. B., & Linden, P. F. (2021). Effects of background rotation on the dynamics of multiphase plumes. *Journal of Fluid Mechanics, 915*, A2. https://doi.org/10.1017/jfm.2020.1181

Gao, Y., Kamenkovich, I., & Perlin, N. (2023). Origins of mesoscale mixed-layer depth variability in the Southern Ocean. *Ocean Science, 19*(3), 615–627. https://doi.org/10.5194/os-19-615-2023

Gaspar, P. (1988). Modeling the seasonal cycle of the upper ocean. *Journal of Physical Oceanography, 18*(2), 161–180. https://doi.org/10.1175/1520-0485(1988)018(0161:MTSCOT)2.0:CO;2

Gaspar, P., Grégoris, Y., & Lefevre, J.-M. (1990). A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station Papa and long-term upper ocean study site. *Journal of Geophysical Research, 95*(C9), 16179–16193. https://doi.org/10.1029/JC095iC09p16179

Gent, P. R. (2011). The Gent–McWilliams parameterization: 20/20 hindsight. *Ocean Modelling, 39*(1–2), 2–9. https://doi.org/10.1016/j.oemod.2010.08.002

Gent, P. R., & McWilliams, J. C. (1990). Isopycnal mixing in ocean circulation models. *Journal of Physical Oceanography, 20*(1), 150–155. https://doi.org/10.1175/1520-0485(1990)020(0150:IMIOCM)2.0:CO;2

Gent, P. R., Willebrand, J., McDougall, T. J., & McWilliams, J. C. (1995). Parameterizing eddy-induced tracer transports in ocean circulation models. *Journal of Physical Oceanography, 25*(4), 463–474. https://doi.org/10.1175/1520-0485(1995)025(0463:PEITTI)2.0:CO;2

Giordani, H., Bourdallé-Badie, R., & Madec, G. (2020). An eddy-diffusivity mass-flux parameterization for modeling oceanic convection. *Journal of Advances in Modeling Earth Systems, 12*(9), e2020MS002078. https://doi.org/10.1029/2020MS002078

Griffies, S. M., Danabasoglu, G., Durack, P. J., Adcroft, A. J., Balaji, V., Böning, C. W., et al. (2016). OMIP contribution to CMIP6: Experimental and diagnosticprotocol for the physical component of the ocean model intercomparison project. *Geoscientific Model Development, 9*(9), 3231–3296. https://doi.org/10.5194/gmd-9-3231-2016

Gutjahr, O., Brüggemann, N., Haak, H., Jungclaus, J. H., Putrasahan, D. A., Lohmann, K., & von Storch, J.-S. (2021). Comparison of ocean vertical mixing schemes in the max Planck Institute Earth system model (MPI-ESM1.2). *Geoscientific Model Development, 14*(5), 2317–2349. https://doi.org/10.5194/gmd-14-2317-2021

Hanjalić, K., & Launder, B. E. (1972). A Reynolds stress model of turbulence and its application to thin shear flows. *Journal of Fluid Mechanics, 52*(4), 609–638. https://doi.org/10.1017/S002211207200268X

Heuzé, C. (2017). North Atlantic deep water formation and AMOC in CMIP5 models. *Ocean Science, 13*(4), 609–622. https://doi.org/10.5194/os-13-609-2017

Holte, J. W., Talley, L. D., Chereskin, T. K., & Sloyan, B. M. (2012). The role of air-sea fluxes in Subantarctic Mode Water formation. *Journal of Geophysical Research, 117*(C3), 2011JC007798. https://doi.org/10.1029/2011JC007798
Johnson, L., Fox-Kemper, B., Li, Q., Pham, H. T., & Sarkar, S. (2023). A finite-time ensemble method for mixed layer model comparison. *Journal of Physical Oceanography*, 53(9), 2211–2230. https://doi.org/10.1175/JPO-D-22-0107.1

Koenigk, T., Fuentes-Franco, R., Meccia, V. L., Gutjahr, O., Jackson, L. C., New, A. L., et al. (2021). Deep mixed ocean volume in the Labrador Sea in HighResMIP models. *Climate Dynamics*, 57(7–8), 1895–1918. https://doi.org/10.1007/s00382-021-05785-x

Kraus, E. B., & Turner, J. S. (1967). A one-dimensional model of the seasonal thermocline II. The general theory and its consequences. *Tellus*, 19(1), 98–106. https://doi.org/10.1111/j.2153-3490.1967.tb01462.x

Kuhlbrodt, T., Griesel, A., Montoya, M., Levermann, A., Hofmann, M., & Rahmstorf, S. (2007). On the driving processes of the Atlantic meridional overturning circulation. *Reviews of Geophysics*, 45(2), RG2001. https://doi.org/10.1029/2004RG000166

Large, W. G., McWilliams, J. C., & Doney, S. C. (1994). Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. *Reviews of Geophysics*, 32(4), 363–403. https://doi.org/10.1029/94RG01872

Lazar, A., Madec, G., & Delecluse, P. (1999). The deep interior downwelling, the veronis effect, and mesoscale tracer transport parameterizations in an ogcm. *Journal of Physical Oceanography*, 29(11), 2945–2961. https://doi.org/10.1175/1520-0485(1999)029(2945:TDIDTV)2.0.CO;2

Legay, A. (2024a). Data for Legay et al. "A framework for assessing ocean mixed layer depth evolution" [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.11236177

Legay, A. (2024b). Legaya/James2024ParameterSpace: 1.0 [Software]. *Zenodo*. https://doi.org/10.5281/zenodo.11381893

Lellouche, J.-M., Chassignet, E., Bourdallé-Badie, R., Garric, G., Melet, A., Drévillon, M., et al. (2021). The Copernicus global 1/12° oceanic and sea ice GLORYS12 reanalysis. *Frontiers in Earth Science*, 9, 698876. https://doi.org/10.3389/feart.2021.698876

Lévy, M., Mémery, L., & Madec, G. (1998). The onset of a bloom after deep winter convection in the northwestern Mediterranean Sea: Mesoscale process study with a primitive equation model. *Journal of Marine Systems*, 16(1–2), 7–21. https://doi.org/10.1016/S0924-7963(97)00097-3

Li, M., & Garrett, C. (1993). Cell merging and the jet/downwelling ratio in Langmuir circulation. *Journal of Marine Research*, 51(4), 737–769. https://doi.org/10.1175/1520-0485(1997)027(0121:MLDDTL)2.0.CO;2

Li, Q., Reichl, B. G., Fox-Kemper, B., Adcroft, A. J., Belcher, S. E., Danabasoglu, G., et al. (2019). Comparing Ocean surface boundary vertical mixing schemes including Langmuir turbulence. *Journal of Advances in Modeling Earth Systems*, 11(11), 3545–3592. https://doi.org/10.1029/2019MS001810

Mack, S. A., & Schoeberlein, H. C. (2004). Richardson number and ocean mixing: Towed chain observations. *Journal of Physical Oceanography*, 34(4), 736–754. https://doi.org/10.1175/1520-0485(2004)034(0736:RNAOMT)2.0.CO;2

Madec, G., Bourdallé-Badie, R., Chanut, J., Clementi, E., Coward, A., Ethé, C., et al. (2022). NEMO ocean engine. *Zenodo Version Number*: v4.2. https://doi.org/10.5281/ZENODO.1464816

Marshall, J., & Schott, F. (1999). Open-ocean convection: Observations, theory, and models. *Reviews of Geophysics*, 37(1), 1–64. https://doi.org/10.1029/98RG02739

Mellor, G. L. (1973). Analytic prediction of the properties of stratified planetary surface 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., & Durbin, P. A. (1975). The structure and dynamics of the ocean surface mixed layer. *Journal of Physical Oceanography*, 5(4), 718–728. https://doi.org/10.1175/1520-0485(1975)005(0718:TSADOT)2.0.CO;2

Mellor, G. L., & Yamada, T. (1974). A hierarchy of turbulence closure models for planetary boundary layers. *Journal of the Atmospheric Sciences*, 31(7), 1791–1806. https://doi.org/10.1175/1520-0469(1974)031(1791:AHOTCM)2.0.CO;2

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

Obukhov, A. M. (1971). Turbulence in an atmosphere with a non-uniform temperature. *Boundary-Layer Meteorology*, 2(1), 7–29. https://doi.org/10.1007/BF00718085

Pollard, R. T., Rhines, P. B., & Thompson, R. O. R. Y. (1973). The deepening of the wind-mixed layer. *Geophysical Fluid Dynamics*, 4(4), 381–404. https://doi.org/10.1080/03091927208236105

Price, J. F. (1979). On the scaling of stress-driven entrainment experiments. *Journal of Fluid Mechanics*, 90(3), 509–529. https://doi.org/10.1017/S0022112079002366

Price, J. F., Weller, R. A., & Pinkel, R. (1986). Diurnal cycling: Observations and models of the upper ocean response to diurnal heating, cooling, and wind mixing. *Journal of Geophysical Research*, 91(C7), 8411–8427. https://doi.org/10.1029/JC091iC07p08411

Reichl, B. G., Adcroft, A., Griffies, S. M., & Hallberg, R. (2022). A potential energy analysis of ocean surface mixed layers. *Journal of Geophysical Research: Oceans*, 127(7), e2021JC018140. https://doi.org/10.1029/2021JC018140

Reichl, B. G., & Hallberg, R. (2018). A simplified energetics based planetary boundary layer (ePBL) approach for ocean climate simulations. *Ocean Modelling*, 132, 112–129. https://doi.org/10.1016/j.ocemod.2018.10.004

Rodi, W. (1987). Examples of calculation methods for flow and mixing in stratified fluids. *Journal of Geophysical Research*, 92(C5), 5305. https://doi.org/10.1029/JC092iC05p05305

Sallée, J.-B., Pellicherero, V., Akhoudas, C., Pauthenet, E., Vignes, L., Schmidtko, S., et al. (2021). Summertime increases in upper-ocean stratification and mixed-layer depth. *Nature*, 591(7851), 592–598. https://doi.org/10.1038/s41586-021-03303-x

Sallée, J. B., Speer, K. G., & Rintoul, S. R. (2010). Zonally asymmetric response of the Southern Ocean mixed-layer depth to the Southern Annular Mode. *Nature Geoscience*, 3(4), 273–279. https://doi.org/10.1038/ngeo812

Schulzweida, U. (2023). Cdo user guide. *Zenodo*. https://doi.org/10.5281/zenodo.10020800

Sérazin, G., Tréguéier, A. M., & De Boyer Montégut, C. (2023). A seasonal climatology of the upper ocean pycnocline. *Frontiers in Marine Science*, 10, 1120112. https://doi.org/10.3389/fmar.2023.1120112

Shy, S. (1995). Mixing dynamics of jet interaction with a sharp density interface. *Experimental Thermal and Fluid Science*, 10(3), 355–369. https://doi.org/10.1016/0894-1777(94)00095-P

Simpson, J. H., & Hunter, J. R. (1974). Fronts in the Irish sea. *Nature*, 250(5465), 404–406. https://doi.org/10.1038/250404a0

Souza, A. N., Wagner, G. L., Ramadhan, A., Allen, B., Churavy, V., Schloss, J., et al. (2020). Uncertainty Quantification of ocean parameterizations: Application to the K-Profile-Parameterization for penetrative convection. *Journal of Advances in Modeling Earth Systems*, 12(12), e2020MS002108. https://doi.org/10.1029/2020MS002108

Speer, K. G., & Marshall, J. (1995). The growth of convective plumes at seafloor hot springs. *Journal of Marine Research*, 53(6), 1025–1057. https://doi.org/10.1357/0022240953212972

Stull, R. B. (Ed.) (1988). *An introduction to boundary layer meteorology*. Springer Netherlands. https://doi.org/10.1007/978-94-009-3027-8

Sverdrup, H. (1953). On conditions for the vernal blooming of phytoplankton. *J Cons Int Explor Mer*, 18(3), 287–295. https://doi.org/10.1093/icesjms/18.3.287

Taylor, J. R., & Ferrari, R. (2011). Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms. *Limnology & Oceanography*, 56(6), 2293–2307. https://doi.org/10.4319/lo.2011.56.6.2293
Taylor, J. R., & Thompson, A. F. (2023). Submesoscale dynamics in the upper ocean. *Annual Review of Fluid Mechanics*, 55(1), 103–127. https://doi.org/10.1146/annurev-fluid-031422-095147

Thangam, S., Abid, R., & Speziale, C. G. (1992). Application of a new K-tau model to near wall turbulent flows. *AIAA Journal*, 30(2), 552–554. https://doi.org/10.2514/3.10952

Thomas, L. N., & Lee, C. M. (2005). Intensification of ocean fronts by down-front winds. *Journal of Physical Oceanography*, 35(6), 1086–1102. https://doi.org/10.1175/JPO2737.1

Treguier, A. M., De Boyer Montégut, C., Bozec, A., Chassignet, E. P., Fox-Kemper, B., McC. Hogg, A., et al. (2023). The mixed-layer depth in the ocean model intercomparison project (OMIP): Impact of resolving mesoscale eddies. *Geoscientific Model Development*, 16(13), 3849–3872. https://doi.org/10.5194/gmd-16-3849-2023

Tsujino, H., Urakawa, S., Nakano, H., Small, R. J., Kim, W. M., Yeager, S. G., et al. (2018). JRA-55 based surface dataset for driving ocean-seaice models (JRA55-do). *Ocean Modelling*, 130, 79–139. https://doi.org/10.1016/j.ocemod.2018.07.002

Turner, J. S. (1973). *Buoyancy effects in fluids*. Cambridge University Press. https://doi.org/10.1017/CBO9780511608827

Turner, J. S. (1986). Turbulent entrainment: The development of the entrainment assumption, and its application to geophysical flows. *Journal of Fluid Mechanics*, 173, 431–471. https://doi.org/10.1017/S0022112086001222

Umlauf, L., & Burchard, H. (2003). A generic length-scale equation for geophysical turbulence models. *Journal of Marine Research*, 61(2), 235–265. https://doi.org/10.1357/002224003322005087

Umlauf, L., & Burchard, H. (2005). Second-order turbulence closure models for geophysical boundary layers. A review of recent work. *Continental Shelf Research*, 25(7–8), 795–827. https://doi.org/10.1016/j.csr.2004.08.004

Ushijima, Y., & Yoshikawa, Y. (2020). Mixed layer deepening due to wind-induced shear-driven turbulence and scaling of the deepening rate in the stratified ocean. *Ocean Dynamics*, 70(4), 505–512. https://doi.org/10.1007/s10236-020-01344-w

Van Der Laan, M. P., Kelly, M., Floors, R., & Peña, A. (2020). Rossby number similarity of an atmospheric RANS model using limited-length-scale turbulence closures extended to unstable stratification. *Wind Energy Science*, 5(1), 355–374. https://doi.org/10.5194/wes-5-355-2020

Van Roekel, L., Adcroft, A. J., Danabasoglu, G., Griffies, S. M., Kauffman, B., Large, W., et al. (2018). The KPP boundary layer scheme for the ocean: Revisiting its formulation and benchmarking one-dimensional simulations relative to LES. *Journal of Advances in Modeling Earth Systems*, 10(11), 2647–2685. https://doi.org/10.1029/2018MS001336

Vreugdenhil, C. A., & Gayen, B. (2021). Ocean convection. *Fluids*, 6(10), 360. https://doi.org/10.3390/fluids6100360

Wagner, G. L., Hillier, A., Constantinou, N. C., Silvestri, S., Souza, A., Burns, K., et al. (2023). *Catke: A turbulent-kinetic-energy-based parameterization for ocean microturbulence with dynamic convective adjustment*. arXiv preprint arXiv:2306.13204.

Weber, J. E. (1983). Steady wind- and wave-induced currents in the open ocean. *Journal of Physical Oceanography*, 13(3), 524–530. https://doi.org/10.1175/1520-0485(1983)013(0524:SWAWIC)2.0.CO;2

Wilcox, D. C. (1988). Reassessment of the scale-determining equation for advanced turbulence models. *AIAA Journal*, 26(11), 1299–1310. https://doi.org/10.2514/3.10041

Willis, G., & Deardorff, J. (1974). A laboratory model of the unstable planetary boundary layer. *Journal of the Atmospheric Sciences*, 31(5), 1297–1307. https://doi.org/10.1175/1520-0469(1974)031(1297:ALMOTU)2.0.CO;2

Wyngaard, J. C. (1973). On the surface-layer turbulence. In *Workshop on micrometeorology* (pp. 101–149).

Zeierman, S., & Wolfshtein, M. (1986). Turbulent time scale for turbulent-flow calculations. *AIAA Journal*, 24(10), 1606–1610. https://doi.org/10.2514/3.9490

Zheng, Z., Harcourt, R. R., & D'Asaro, E. A. (2021). Evaluating Monin–Obukhov scaling in the unstable oceanic surface layer. *Journal of Physical Oceanography*, 51(3), 911–930. https://doi.org/10.1175/JPO-D-20-0201.1

Zhu, Y., Zhang, R.-H., & Sun, J. (2020). North Pacific upper-ocean cold temperature biases in CMIP6 simulations and the role of regional vertical mixing. *Journal of Climate*, 33(17), 7523–7538. https://doi.org/10.1175/JCLI-D-19-0654.1