# A Generalized Stochastic Formulation of the
Ekman–Stokes Model with Statistical Analyses

Long Li, Etienne Mémin, Bertrand Chapron

► To cite this version:

Long Li, Etienne Mémin, Bertrand Chapron. A Generalized Stochastic Formulation of the Ekman–Stokes Model with Statistical Analyses. Journal of Physical Oceanography, 2025, 55 (9), pp.1389-1407. 10.1175/JPO-D-24-0116.1 . hal-05327658

HAL Id: hal-05327658

https://inria.hal.science/hal-05327658v1

Submitted on 3 Dec 2025

**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
A Generalized Stochastic Formulation of the Ekman–Stokes Model with
Statistical Analyses

LONG LI, a,b ETIENNE MÉMIN, a,b AND BERTRAND CHAPRON a,c

a Odyssey Team, Centre Inria de l'Université de Rennes, Rennes, France

b IRMAR, Université de Rennes, Rennes, France

c Ifremer, Univ. Brest, CNRS, IRD, Laboratoire d'Océanographie Physique et Spatiale, IUEM, Plouzané, France

(Manuscript received 26 June 2024, in final form 9 April 2025, accepted 12 May 2025)

**ABSTRACT:** To describe the upper-ocean Ekman boundary layer, a novel stochastic model is formulated to better capture the interplay between random components associated with the numerically unresolved physical process induced by wind, waves, and currents. Accounting for the uncertainty of unresolved velocity fluctuations, resulting Ekman velocities display increased variability, higher kinetic energy, and more frequent extreme events than in an Ekman-Stokes model with time-correlated additive Gaussian forcing. Near the surface, stronger correlations and skewed probability distributions are also revealed. Sensitivity analyses highlight impacts of transient winds and surface waves, which deepen circulation and increase the dispersion of realizations in a balanced way according to error representations. Energy and transport are maximized when the mean directions of wind velocity and surface waves coincide. This model offers a new approach to describing upper-ocean dynamics, providing improved insights for anticipating vertical fluxes, and enhancing predictive capabilities.

SIGNIFICANCE STATEMENT: The upper-ocean boundary layer plays a critical role in regulating air-sea interactions, vertical mixing, and climate variability. Many existing boundary layer models employ deterministic turbulence closures that do not account for the inherently random nature of unresolved small-scale processes. This study presents a novel stochastic framework that explicitly incorporates the random effects of wind forcing, surface waves, and turbulent transport, while maintaining key physical balances. The results demonstrate enhanced variability and a higher frequency of extreme events, emphasizing the relevance of stochasticity in upper-ocean dynamics. This approach provides a physically consistent representation of boundary layer processes, with potential benefits for improving ocean modeling and forecasting of an ensemble of realizations at coarse resolution.

KEYWORDS: Boundary layer; Statistics; Uncertainty; Dynamical system model; Stochastic models

**1. Introduction**

The seminal work of Ekman (1905) provides an archetypi-
cal exact solution of the simplified Navier-Stokes equations.
It is fundamental to the theory of ocean circulation and attrac-
tive for the theoretical analysis. Key assumptions are to con-
sider the upper-ocean boundary layer flow steady, linear, of
uniform density and viscosity, and driven solely by a surface
wind stress. A resulting balance between Coriolis and turbu-
lent drag forces shall then lead to the gradual rotation and de-
cay of ageostrophic velocity with depth, presenting a spiral
vertical structure. Among the anticipated limitations of this
model, the definition of turbulent viscosity—used to parame-
terize, exclusively through a diffusion operator, the effect of
unresolved velocity components due to turbulent motion or
small eddies—has been called into question. Beyond possibly
contributing to nonuniform vertical mixing (Large et al.
1994), the unresolved fluctuating motions redistribute and en-
hance energy in a manner that statistically preserves the

fluctuation-dissipation balance. This latter contribution is, in
general, not taken into account in many existing models that
rely on deterministic turbulence closures.

Significant modifications of the Ekman spiral may also re-
sult from near-surface currents induced by surface gravity
waves (Huang 1979; Jenkins 1986; Xu and Bowen 1994;
McWilliams et al. 1997, 2012), not all necessarily aligned with
the wind direction. Buoyancy effects can also play a crucial
role in the boundary layer dynamics, particularly at fronts
(Price and Sundermeyer 1999; McWilliams et al. 2009; Gula
et al. 2014; McWilliams et al. 2015). Although incorporating
these additional effects complicates analytical solutions, ad-
vanced mathematical tools (Lewis and Belcher 2004; Wenegrat
and McPhaden 2016b,a; Higgins et al. 2020) and accurate nu-
merical methods offer viable approximations.

Generalized Ekman-Stokes models have been proposed to
account for surface wave effects through the Stokes drift
(McWilliams et al. 1997) and the variability of wind stress
(McWilliams and Huckle 2006). In homogeneous mixing has
also been considered through manageable forms such as the
nonlocal K-profile parameterization (KPP) (McWilliams and
Huckle 2006; Wenegrat and McPhaden 2016a). More gener-
ally, the Ekman balance can be explicitly reformulated as a
stochastic differential equation by introducing additive ran-
dom terms that are either time correlated (McWilliams et al.

<footnote>Denotes content that is immediately available upon publica-
tion as open access.</footnote>

<footnote>Corresponding author: Long Li, long.li@inria.fr</footnote>
1997; McWilliams and Huckle 2006) or decorrelated (Chu 2008), to represent wind forcing. However, in all the above developments, the statistically consistent impact of small-scale unresolved processes—particularly their nondiffusive contri- butions to large-scale dynamics—remains not accounted for, and the fluctuation–dissipation imbalance still persists. The generalized models effectively capture only the average dissi- pation effects associated with unresolved motions as well as their induced advective correction (e.g., the Stokes drift) but do not explicitly describe the actions of small-scale processes. The complexities introduced by non-Gaussian perturbations, such as extreme and transient events—e.g., wave groupiness, breaking events, and intermittency—are potentially completely neglected.

Extremely challenging, these effects may be difficult, not impossible, to tackle within a deterministic one-dimensional (1D) vertical model. In contrast, following a stochastic framework, the specification of these interactions can be considered to enable their explicit integration into the model. To our knowledge, none of the proposed stochastic models mentioned previously include a stochastic parameterization of the small-scale fluctuations associated both with wind forcing, surface waves, and turbulence. To be relevant, such a model should adhere to the scaling laws and derivations that yield the 1D velocity profile of the Ekman boundary layer, while incorporating at its core a random representation of the three-dimensional (3D) velocity fluctuations within a framework guaranteeing a fluctuation–dissipation balance.

In this study, we propose to rely on a recently developed stochastic formalism, which has notably been applied to extend the classical Craik–Leibovich model of Langmuir turbulence (Bauer et al. 2020). Originally proposed in Mémin (2014), and referred to as modeling under location uncertainty (LU), such a formalism has demonstrated versatility in the derivation of a wide range of models, including reduced-order models (Chapron et al. 2018), quasigeostrophic models (Li et al. 2023a), wave models (Mémin et al. 2024), and stochastic primitive hydrostatic models (Tucciarone et al. 2024). LU-based models strictly respect the fluctuation–dissipation balance and have been shown to effectively reproduce long-term statistical properties (Li et al. 2023a) and capture large-scale flow structuration (Bauer et al. 2020).

Owing to its foundation in the stochastic interpretation of the Reynolds transport theorem, the LU framework preserves the theoretical energy conservation properties of deterministic systems (e.g., energy or tracer moments) (Resseguier et al. 2017; Brecht et al. 2021; Li et al. 2023a). For the Navier–Stokes equations, it also retains the same asymptotic convergence properties (Debussche et al. 2023), such as uniqueness in 2D and the existence of weak probabilistic solutions in 3D. These results arise from a strict balance between energy input from stochastic forcing and system dissipation—that is, LU models directly respect the fluctuation–dissipation theorem (Bauer et al. 2020; Resseguier et al. 2017). Moreover, LU models converge to the deterministic dynamics as noise vanishes (Debussche et al. 2023). This convergence is nontrivial: Even weak, ad hoc non-Gaussian noise can significantly alter the distribution of velocity increments (Chapron et al. 2018).

To describe upper-ocean balance, the LU framework allows us to consistently integrate the combined effects of random winds, surface waves, and turbulent mixing. We build on this ability to represent not only the average impact of these processes but also their impulsive, potentially non-Gaussian, and correlated behaviors. This approach supports exploration of the joint influence of surface waves, wind stress, and turbulent mixing on the Ekman model. Stratification effects are omitted here and will be addressed in future work.

To represent small-scale unresolved components, random fields can be constructed such that their covariances are consistent with established vertical mixing closures, such as the KPP viscosity and the Stokes drift. This approach enables a direct connection to existing generalized Ekman models (Wenegrat and McPhaden 2016a,b) and Ekman–Stokes formulations (McWilliams et al. 1997; Higgins et al. 2020). As already proposed by McWilliams and Huckle (2006), wind stress may be modeled as an external random forcing via an Ornstein–Uhlenbeck process.

In this paper, by appropriately scaling the generalized stochastic Craik–Leibovich equation outlined in Bauer et al. (2020), both steady and time-dependent stochastic Ekman–Stokes models are derived. An approximate analytical solution to the steady model is presented using a well-established approach (Wenegrat and McPhaden 2016b). The time-dependent model, formulated in one dimension, is investigated through large-ensemble simulations to quantify uncertainties and system responses to various random inputs (wind, waves, and turbulence), with particular emphasis on the occurrence of extreme events. The statistical responses of the proposed stochastic model are then compared to those of a deterministic Ekman–Stokes model with random wind forcing, followed by a sensitivity analysis of the stochastic model’s response to varying wind and surface wave parameters.

The rest of this paper is structured as follows: Section 2 presents the proposed stochastic Ekman–Stokes models with detailed derivations and describes the uncertainty parameterization methods. Section 3 discusses the numerical results, including the statistical diagnosis and sensitivity analysis of the random models. Finally, in section 4, we draw conclusions and provide an outlook for future research endeavors.

## 2. Stochastic formulation

In this section, we first briefly review the generalized stochastic Craik–Leibovich equations as derived in our previous work (Bauer et al. 2020). Stochastic Ekman–Stokes models are then obtained by appropriately scaling these nonlinear equations. Subsequently, we describe the representation of uncertainties induced by established physical parameterizations.

### a. Ekman-layer scaling in stochastic Craik–Leibovich equations

A general open question in fluid dynamics concerns how small-scale processes influence the large-scale structure of fluid flows. In oceanic and atmospheric contexts, wave–current interactions—particularly Langmuir circulation—illustrate this effect. These interactions, modeled by the Craik–Leibovich system, describe how the wave-driven Stokes drift generates a
"vortex force," resulting in streak formations within the flow. In our previous work (Bauer et al. 2020), we showed that this effect can be generalized as a consequence of the statistical inhomogeneity of the small-scale flow component through a stochastic representation of the flow.

The adopted stochastic framework, termed location uncertainty (Mémin 2014), arises from a decomposition of the Lagrangian fluid flow into a resolved, time-smooth flow component and an unresolved, highly oscillatory noise term. This noise term,

denoted by $\sigma dB_t$, is represented as a centered, time-uncorrelated but spatially correlated Gaussian process, where $B$ is a cylindrical Brownian motion (Da Prato and Zabczyk 2014) defined on an infinite-dimensional Hilbert space and $\sigma$ is a random correlation process (Debussche et al. 2023; Li et al. 2023a,b) defined within a space of Hilbert-Schmidt operators and mapping two spaces of functions with bounded variance. Within this framework, the generalized stochastic Craik-Leibovich momentum equations (Bauer et al. 2020) are expressed as follows:

$$ 
\begin{aligned}
dv + & \left[ \mathbf{v} \cdot \nabla \mathbf{v} + f \hat{\mathbf{z}} \times \mathbf{v} + \underbrace{\mathbf{v}_s \times \mathbf{v}_s}_{\text{Coriolis-Stokes force}} + \underbrace{\mathbf{v}_s \times \mathbf{v}_s}_{\text{Vortex force}} - b \hat{\mathbf{z}} + \nabla \left( \underbrace{p + \mathbf{v} \cdot \mathbf{v}_s}_{\text{Effective pressure}} \right) \right] dt \\
& = \left\{ \begin{array}{c} \underbrace{\nabla \cdot [\mathbf{a}\nabla(\mathbf{v} + \mathbf{v}_s)]}_{\text{Subgrid-scale diffusion}} + \underbrace{\mathbf{v}_s \cdot \nabla \mathbf{v}}_{\text{Itô-Stokes advection}} - \underbrace{\mathbf{v}_s \times \mathbf{v}}_{\text{Itô-Stokes force}} \\ \left\{ \begin{array}{c} \underbrace{\nabla \cdot [\mathbf{a}\nabla(\mathbf{v} + \mathbf{v}_s)]}_{\text{Subgrid-scale diffusion}} + \underbrace{\mathbf{v}_s \cdot \nabla \mathbf{v}}_{\text{Itô-Stokes advection}} - \underbrace{\mathbf{v}_s \times \mathbf{v}}_{\text{Itô-Stokes force}} \end{array} \right\} dt - \left[ \begin{array}{c} \underbrace{\sigma \cdot \nabla(\mathbf{v} + \mathbf{v}_s)}_{\text{Subgrid-scale advection}} + f \hat{\mathbf{z}} \times \sigma \end{array} \right] dB_t - \nabla dP_t \end{aligned} 
\quad (2.1a) $$

$$ \nabla \cdot \mathbf{v} = 0, \nabla \cdot \sigma dB_t = 0, \quad (2.1b) $$

$$ \mathbf{a} := \frac{1}{2} \boldsymbol{\sigma} \boldsymbol{\sigma}^{\mathrm{T}}, \mathbf{v}_s := \nabla \cdot \mathbf{a}. \quad (2.1c) $$

Here, $\mathbf{v} = (u, v, w)^{\mathrm{T}}$ denotes the resolved effective velocity, which incorporates the statistical effects of unresolved noise, and $dv$ represents its time evolution at a fixed point. The three-dimensional gradient operator is $\nabla = (\partial_x, \partial_y, \partial_z)^{\mathrm{T}}$, $f$ represents the Coriolis frequency, and $\hat{\mathbf{z}} = (0, 0, 1)^{\mathrm{T}}$ represents the vertical unit vector. The vorticity, $\boldsymbol{\omega} = \nabla \times \mathbf{v}$, is defined as the curl of the velocity. The buoyancy variable $b = gp'/\rho_o$ is introduced under the Boussinesq approximation, where $g$ is the gravitational constant, $\rho_o$ is the reference density of water, and $\rho'$ is the (small) density fluctuation. The resolved dynamic pressure $p$ is scaled by $\rho_o$.

Two statistical properties [Eq. (2.1c)] of the noise play a critical role in the preceding stochastic system. The random, symmetric, nonnegative diffusion tensor $\mathbf{a}$ (in units of $m^2 s^{-1}$), referred to as Itô diffusion tensor, is defined by the quadratic covariation of the noise process (Bauer et al. 2020). It can be interpreted as the noise's conditional variance with respect to the correlation operator, scaled by the time step: $\mathbf{a} = (1/2)\mathbb{E}[\boldsymbol{\sigma}dB_t(\boldsymbol{\sigma}dB_t)^{\mathrm{T}}|\boldsymbol{\sigma}]/dt = (1/2)\boldsymbol{\sigma}\boldsymbol{\sigma}^{\mathrm{T}}$, where $\mathbb{E}$ represents the expectation value. The Itô-Stokes drift velocity, $\mathbf{v}_s = \nabla \cdot \mathbf{a}$, is a statistical effect that arises from the inhomogeneity of the unresolved random field. It corresponds to the turbophoresis phenomenon, wherein inertial particles migrate toward regions of lower turbulent diffusivity (Reeks 1983). When the noise is homogeneous—meaning that $\boldsymbol{\sigma}$ acts through a spatial convolution operator and can be diagonalized in Fourier modes (Mémin 2014; Mémin et al. 2024)—the diffusion tensor $\mathbf{a}$ becomes a constant matrix at each point, causing the statistical drift $\mathbf{v}_s$ to vanish. Notably, both $\mathbf{a}$ and $\mathbf{v}_s$ arise naturally within the framework of Itô's stochastic calculus (Bauer et al. 2020). This formulation assumes a quasi-stationary $\mathbf{v}_s$ for the derivation of Eq. (2.1a).

We now interpret several key terms in this stochastic system. The first line reflects the structure of the deterministic Craik-Leibovich (CL) momentum equations (Craik and Leibovich 1976; Leibovich 1980; McWilliams et al. 1997), with the last three terms corresponding to the Coriolis-Stokes force, the CL vortex force, and the Stokes-corrected pressure, respectively. In the stochastic framework, however, the vortex force represents the statistical influence of noise inhomogeneity on large-scale flow dynamics. As a result, this momentum equation provides a generalized stochastic extension of the CL system, in which the statistical drift $\mathbf{v}_s$ replaces the Stokes drift induced by surface wave motion. This distinction broadens the interpretation of $\mathbf{v}_s$, termed Itô-Stokes drift, as it originates not only from wave motion but also from inhomogeneous, unresolved small-scale fluctuations.

The right-hand side of Eq. (2.1a) introduces additional contributions specific to the stochastic formulation. The subgrid-scale (SGS) advection term, $\boldsymbol{\sigma} \cdot \nabla(\mathbf{v} + \mathbf{v}_s)dB_t$, and the SGS diffusion term, $\nabla \cdot [\mathbf{a}\nabla(\mathbf{v} + \mathbf{v}_s)]$, represent the stochastic transport of the resolved momentum ($\mathbf{v} + \mathbf{v}_s$) due to unresolved fluid motions. Within these SGS terms, coupling effects between the Itô-Stokes drift and unresolved turbulence emerge, expressed through $(\boldsymbol{\sigma} \cdot \nabla\mathbf{v}_s)dB_t$ and $\nabla \cdot (\mathbf{a}\nabla\mathbf{v}_s)$. When the Itô-Stokes drift corresponds to a true Stokes drift, as discussed in subsequent sections, these terms account for surface-wave-induced mixing effects. Notably, the diffusive term $\nabla \cdot (\mathbf{a}\nabla\mathbf{v}_s)$ has already been examined in deterministic settings (McWilliams and Sullivan 2000).

Additionally, two further contributions from the Itô-Stokes drift emerge: a corrective advection term, $\mathbf{v}_s \cdot \nabla\mathbf{v}$, and the so-called Itô-Stokes force, $\boldsymbol{\omega}_s \times \mathbf{v}$. The latter term reflects the interaction between the effective velocity and the vorticity of the Itô-Stokes drift, which can be combined with the Coriolis force to form a corrective term. The pressure noise $dP_t$ acts as a Lagrange multiplier enforcing the incompressibility of the unresolved flow through the divergence-free constraints in Eq. (2.1b), which also ensure the energy conservation of the derived stochastic system (Brecht et al. 2021; Li et al. 2023a).
In this study, we consider a simplified stochastic CL system tailored for the Ekman boundary layer, incorporating the effects of surface waves and uncertainties arising from small-scale dynamics. The derivation of the stochastic Ekman-Stokes model is based on scaling the nonlinear stochastic momentum equations [Eq. (2.1a)] and constraining the uncertainties using established physical parameterizations. To isolate the effects of surface waves, we neglect stratification within the boundary layer, assuming a constant buoyancy (Vallis 2017). Without loss of generality, we set $b = 0$ in the subsequent analysis. Key aspects of the derivation are outlined here, with detailed steps provided in appendix A. Most variables follow the same scaling principles as in the deterministic case (Pedlosky 1990). In particular, the horizontal and vertical components of the noise are scaled as

$$ \sigma_x dB_t = \sqrt{2TA_h^{1/2}} \hat{\sigma}_x d\hat{B}_t, \quad \sigma_z dB_t = \sqrt{2TA_v^{1/2}} \hat{\sigma}_z d\hat{B}_t, \quad (2.2a) $$

where $\sigma dB_t = (\sigma_x dB_t, \sigma_z dB_t)^T$ with $\sigma_x = (\sigma_x, \sigma_z)^T$, $T = L/U$ represents the advection time scale, and $L$ and $U$ are the characteristic length and velocity scales, respectively. Here, $A_h$ and $A_v$ denote the horizontal and vertical eddy viscosity coefficients, while the hatted variables represent nondimensional quantities. As a result, the components of the corresponding Itô diffusion tensor are scaled as

$$ \mathbf{a}_{xx} = A_h \hat{\mathbf{a}}_{xx}, \quad a_{zz} = A_v \hat{a}_{zz}, \quad \mathbf{a}_{xz} = A_h^{1/2} A_z^{1/2} \hat{\mathbf{a}}_{xz}. \quad (2.2b) $$

A nondimensional version of the stochastic momentum equations [Eq. (A5)] can be derived with these scalings, along with several classical dimensionless numbers: the Rossby number $Ro = U/(f_0 L)$, the aspect ratio of the horizontal and vertical length scales $\gamma = H/L$, and the horizontal and vertical Ekman numbers $E_h = A_h'(f_0 L^2)$ and $E_v = A_v'(f_0 H^2)$, where $f_0$ and $H$ represent the Coriolis frequency and vertical scales, respectively. To further simplify and derive the equations of motion for the Ekman boundary layer, we adopt the following assumptions regarding the dimensionless numbers:

$$ Ro \ll 1, \quad \gamma \ll 1, \quad E_v \sim \hat{\delta}_e^2, \quad E_h \sim Ro^2, \quad (2.3) $$

where $\hat{\delta}_e$ denotes the nondimensional boundary layer depth. We recall that a small Rossby number suppresses the nonlinear advection terms, while a small aspect ratio ensures hydrostatic balance in the vertical. The third assumption sets the vertical friction to be as large as the Coriolis force in the boundary layer. The last assumption emerges from the ratio of horizontal friction force to inertial acceleration $E_h/Ro = 1/Re$, where $Re = UL/A_h$ is the Reynolds number of the interior flow. In most cases of geophysical interest, the Reynolds number is quite large (Pedlosky 1990). In particular, we choose $1/Re \sim Ro$ here.

Using the above assumptions, the simplified nondimensional momentum equations are expressed as (appendix A)

$$ \begin{aligned} \mathrm{Rod}\hat{\mathbf{u}} = & \{-\hat{f}(\hat{\mathbf{u}} + \partial_{\hat{z}}\hat{\mathbf{a}}_{xz})^{\perp} - \hat{\nabla}\hat{p} + \partial_{\hat{z}}[\hat{a}_{zz}\partial_{\hat{z}}(\hat{\mathbf{u}} + \partial_{\hat{z}}\hat{\mathbf{a}}_{xz})] \\ & + \mathcal{O}(Ro)\}d\hat{t} - \sqrt{2}\mathrm{Ro}^{1/2}\{[\hat{\sigma}_z\partial_{\hat{z}}(\hat{\mathbf{u}} + \partial_{\hat{z}}\hat{\mathbf{a}}_{xz}) \\ & + \hat{f}\hat{\sigma}_x^{\perp}]d\hat{B}_t + \hat{\nabla}d\hat{P}_t\} + \mathcal{O}(Ro)d\hat{B}_t, \end{aligned} \quad (2.4a) $$

$$ \gamma^2 \mathrm{Rod}\hat{w} = \left[ -\frac{1}{\hat{\delta}_e^{1/2}} \partial_{\hat{z}} \hat{p} + \mathcal{O}(\gamma^2) \right] d\hat{t} - \sqrt{2} \frac{\mathrm{Ro}^{1/2}}{\hat{\delta}_e^{1/2}} \partial_{\hat{z}} d\hat{P}_t \\ + \mathcal{O}(\gamma^2) d\hat{B}_t. \quad (2.4b) $$

From this point forward, $\nabla = (\partial_{x_i} \partial_{y_i})^T$ denotes the horizontal gradient (with $\hat{\nabla}$ denoting gradients with respect to adimensional variables and $\partial_{\hat{z}}$ denoting the partial derivative with respect to vertical direction $\hat{z} = \hat{z}/(\hat{\delta}_e)$, $\mathbf{u} = (u, v)^T$ represents the horizontal velocity, and $\mathbf{u}^{\perp} = (-v, u)^T$ stands for its perpendicular counterpart. It is important to note that the nondimensional time variable $\hat{t}$ is not necessarily of $O(1)$. In the following sections, we first derive a steady Ekman-Stokes model, followed by a time-dependent stochastic model, each corresponding to different time scales.

b. Steady Ekman-Stokes model with random coefficients

To derive a steady solution, we consider an intermediate time scale $\hat{t} = O(1)$. The temporal variation of horizontal momentum is neglected; Eq. (2.4a) is thus split into a prognostic equation of $O(1)$ in terms of “d$t$” and a diagnostic equation of $O(Ro^{1/2})$ in terms of “d$\hat{B}_t$.” This splitting, based on physical scaling arguments, is also mathematically justified by the canonical decomposition of a semimartingale (Le Gall 2016). Similarly, the vertical momentum equation [Eq. (2.4b)] reduces to two hydrostatic balances for the resolved and unresolved scales. These are

$$ \hat{f}(\hat{\mathbf{u}} + \hat{\mathbf{u}}_s)^{\perp} = -\hat{p} + \partial_{\hat{z}}[\hat{a}_{zz}\partial_{\hat{z}}(\hat{\mathbf{u}} + \hat{\mathbf{u}}_s)], \quad \hat{\mathbf{u}}_s = \partial_{\hat{z}}\hat{\mathbf{a}}_{xz}, \quad (2.5a) $$

$$ [\hat{f}\hat{\sigma}_x^{\perp} - \hat{\sigma}_z\partial_{\hat{z}}(\hat{\mathbf{u}} + \hat{\mathbf{u}}_s)]d\hat{B}_t = -\hat{\nabla}d\hat{P}_t, \quad (2.5b) $$

$$ \partial_{\hat{z}} \hat{p} = 0, \quad \partial_{\hat{z}} d\hat{P}_t = 0. \quad (2.5c) $$

Following Vallis (2017), assuming $\hat{p}$ independent of $z$ facilitates the separation between pressure-driven interior geostrophic motions and boundary layer ageostrophic motions. Let $\hat{\mathbf{u}} = \hat{\mathbf{u}}_g + \hat{\mathbf{u}}_e$ and $\hat{p} = \hat{p}_g + \hat{p}_e$, where the Ekman layer corrections, denoted with a subscript $e$, are negligible away from the boundary layer. Combining this with the hydrostatic balance, we have $\hat{p}_e = 0$ everywhere, implying there is no boundary layer in the pressure field. Restoring the dimensional form, the ageostrophic velocity in the Ekman layer satisfies

$$ f(\mathbf{u}_e + \mathbf{u}_s)^{\perp} = \partial_z[a_{zz}\partial_z(\mathbf{u}_e + \mathbf{u}_s)], \quad \mathbf{u}_s = \partial_z\mathbf{a}_{xz}. \quad (2.6) $$

The resulting horizontal Itô-Stokes drift in the Ekman boundary layer is given by the vertical derivative of the co-variation process $(\mathbf{a}_{xz})$ between the horizontal and vertical components of the unresolved motions. The previous vectorial equations can be rewritten as a scalar equation for the complex velocity $\mathbf{u} = u + iv$, where $i$ represents the imaginary unit. Together with classical boundary conditions for the momentum flux, the steady Ekman-Stokes model is expressed as

$$ \partial_z(a_{zz}\partial_z\mathbf{u}_v) - if\mathbf{u}_e = if\mathbf{u}_s - \partial_z(a_{zz}\partial_z\mathbf{u}_s), \quad (2.7a) $$

$$ \rho_o a_{zz} \partial_z \mathbf{u}_e |_{z=0} = \tau_w, \quad \rho_o a_{zz} \partial_z \mathbf{u}_e |_{z=-H} = 0, \quad (2.7b) $$
where $\tau_w$ denotes the surface wind stress. In the general case, this elliptic partial differential equation (PDE) involves random coefficients such as the diffusion tensor $\mathbf{a}$, the Itô–Stokes drift $\mathbf{u}_s$, and the wind stress $\tau_w$. It is important to note that the derived formulation (2.7) incorporates the vertical diffusive mixing effect of the Itô–Stokes drift and represents a generalization of the deterministic steady Ekman model with Coriolis–Stokes force.

Following the method proposed by Wenegrat and McPhaden (2016b), an approximate analytical solution can be inferred for Eq. (2.7). To this end, we first vertically differentiate Eq. (2.7a) and then multiply by $\rho_o a_{zz}$ to form an equation for the Ekman current stress $\tau_e = \rho_o a_{zz} \partial_z \mathbf{u}_e$, which reads

$$ (a_{zz} \partial_{zz}^2 - if) \tau_e = -(a_{zz} \partial_{zz}^2 - if) \tau_s, \quad (2.8a) $$

$$ \tau_e(0) = \tau_w, \quad \tau_e(-H) = 0. \quad (2.8b) $$

Similarly, $\tau_s = \rho_o a_{zz} \partial_z \mathbf{u}_s$ represents the derived Itô–Stokes stress. Wenegrat and McPhaden (2016b) propose solving their generalized Ekman model, which shares the same left-hand side differential operator with Eq. (2.8a) but with different right-hand side forces, by initially approximating a solution to the homogeneous equation using the Wentzel–Kramer–Brillouin (WKB) method (Bender and Orszag 1978). Then, the authors solve for the inhomogeneous solution using the method of variation of parameters. Following their detailed derivation, we apply the same approach to solve our boundary value problem (2.8). The complete solution for the Ekman stress, accounting for an additional inhomogeneous forcing, is given by

$$ \tau_e(z) = \tau_w \left[ \frac{a_{zz}(z)}{a_{zz}(0)} \right]^{1/4} \frac{\sinh[\theta(z)]}{\sinh[\theta(0)]} + \int_{-H}^{0} G(z, \zeta) [(if a_{zz}^{-1} - \partial_{zz}^2)\tau_s](\zeta) d\zeta, \quad (2.9a) $$

$$ \theta(z) = \sqrt{if} \int_{-H}^{z} a_{zz}^{-1/2}(\zeta) d\zeta. \quad (2.9b) $$

Here, $G$ denotes a symmetric Green function defined by

$$ G(z, \zeta) = \frac{a_{zz}^{1/4}(z) a_{zz}^{1/4}(\zeta)}{\sqrt{if} \sinh[\theta(0)]} g(z, \zeta), \\ g(z, \zeta) = \begin{cases} \sinh[\theta(z)] \sinh[\theta(\zeta) - \theta(0)] & \text{if } \zeta > z, \\ \sinh[\theta(\zeta)] \sinh[\theta(z) - \theta(0)] & \text{if } \zeta < z. \end{cases} \quad (2.9c) $$

The accuracy of the WKB approximate solution has been further discussed in Wenegrat and McPhaden (2016b), showing generally small errors for small vertical Ekman number $E_v$. Once the Ekman stress is found, the Ekman current solution satisfying Eq. (2.7) can be recovered by

$$ \mathbf{u}_e = -\frac{i}{\rho_o f} \partial_z (\tau_e + \tau_s) - \mathbf{u}_s. \quad (2.10) $$

Integrating this relation vertically, a modified relation for the Ekman transport, $\mathbf{T}_e = \int_{-H}^{0} \mathbf{u}_e(z) dz$, is obtained:

$$ \mathbf{T}_e = -\frac{i}{\rho_0 f} (\tau_w + \tau_s^0) - \mathbf{T}_s. \quad (2.11) $$

This shows that the Ekman transport $\mathbf{T}_e$ is modified by both the surface Itô–Stokes stress $\tau_s^0 := \tau_s(0)$ and the Itô–Stokes transport $\mathbf{T}_s = \int_{-H}^{0} \mathbf{u}_s(z) dz$. Nevertheless, the Ekman–Stokes transport $\mathbf{T}_e + \mathbf{T}_s$ is rotated 90° to the left of the effective surface stress $\tau_w + \tau_s^0$, hence satisfying the deterministic integral relation. Moreover, the vertical velocity induced by the Ekman layer, usually referred to as the Ekman pumping, can be derived by integrating the continuity equation [Eq. (2.1b)] for the ageostrophic component, namely,

$$ w_e(-H) = \frac{1}{\rho_o f} \nabla \times (\tau_w + \tau_s^0) - \nabla \cdot \mathbf{T}_s. \quad (2.12) $$

Accordingly, the Ekman pumping is now modified by both the curl of the surface Itô–Stokes stress and the horizontal divergence of the Itô–Stokes transport. Through the Itô–Stokes drift, it incorporates the statistical effect of the inhomogeneity of the small-scale velocity component. It captures wavy motions, as well as other physical processes affecting small-scale velocity fluctuations, including the effect of vertical mixing.

### c. Time-dependent stochastic Ekman–Stokes model

In this case, a fast time scale is considered by applying a change of time coordinates (Crowe and Taylor 2018), such that $\hat{\tau} = \hat{f}/\text{Ro}$, and hence $\hat{B}_\tau = \hat{B}_t/\text{Ro}^{1/2}$. Therefore, Eqs. (2.4a) and (2.4b) reduce to

$$ d\hat{\mathbf{u}} = \{-\hat{f}(\hat{\mathbf{u}} + \hat{\mathbf{u}}_s)^\perp - \hat{\nabla}\hat{p} + \partial_\zeta[\hat{a}_{zz}\partial_\zeta(\hat{\mathbf{u}} + \hat{\mathbf{u}}_s)]\}d\hat{\tau} - \sqrt{2}[\hat{f}\hat{\sigma}_x^\perp + \hat{\sigma}_z\partial_\zeta(\hat{\mathbf{u}} + \hat{\mathbf{u}}_s)]d\hat{B}_\tau - \sqrt{2}\hat{\nabla}d\hat{P}_\tau, \quad (2.13a) $$

$$ \partial_\zeta \hat{p} = 0, \quad \partial_\zeta d\hat{P}_\tau = 0. \quad (2.13b) $$

Unlike the steady case, the contributions of the $d\hat{t}$ and $d\hat{B}_t$ terms on the momentum evolution can no longer be separated. Both terms have exactly the same scales. By separating the pressure-driven geostrophic motions, the wind-driven ageostrophic motions in the Ekman layer satisfy the following dimensional stochastic PDE (SPDE) (in complex notations):

$$ d\mathbf{u}_e = \{-if(\mathbf{u}_e + \mathbf{u}_s) + \partial_z[a_{zz}\partial_z(\mathbf{u}_e + \mathbf{u}_s)]\}dt - [if\boldsymbol{\sigma}_x + \sigma_z\partial_z(\mathbf{u}_e + \mathbf{u}_s)]dB_t, \quad (2.14a) $$

$$ \rho_o a_{zz} \partial_z \mathbf{u}_e |_{z=0} = \tau_w, \quad \rho_o a_{zz} \partial_z \mathbf{u}_e |_{z=-H} = 0. \quad (2.14b) $$

This stochastic formulation extends the deterministic time-dependent Ekman–Stokes model (McWilliams et al. 1997), incorporating the vertical unresolved advective process of the Ekman current, the random unresolved Coriolis force, and both the diffusive and advective vertical mixing effects of the Itô–Stokes drift. This system has the advantage of not only incorporating large-scale vertical mixing terms but also accounting for unresolved vertical advection contributing to it. By this way, one can transition from an analysis in mean to an analysis of higher-order statistical moments. The role of the
Stokes drift, when interpreted as a special case of the Itô–Stokes drift (as will be the case in the following), is considered both in the mixing term and in the fluctuation term. Its impli- cations on the solution will be discussed in detail in the nu- merical results section.

Finding analytical expression of strong (pathwise) solutions to an SPDE with multiplicative noises is extremely challeng- ing. Instead, in section 3, the statistical moments for solutions of the SPDE (2.14), incorporating other random parameters is investigated through numerical simulations based on the Monte Carlo method (Caflisch 1998).

Here, we exemplify important integral properties of the conditional mean solution. Assuming that the random corre- lation operator σ and the cylindrical Brownian motion B are independent, and taking the expectation of Eq. (2.14a) with respect to the conditional probability distribution given σ , we deduce the following PDE for E $^{σ}$[ u$_{e}$ ] := E [ u$_{e}$ | σ ]:

$$
\partial_t \mathbb{E}^\sigma[\mathbf{u}_e] = -if(\mathbb{E}^\sigma[\mathbf{u}_e] + \mathbf{u}_s) + \partial_z[a_{zz}\partial_z(\mathbb{E}^\sigma[\mathbf{u}_e] + \mathbf{u}_s)]. \quad (2.15)
$$

This resulting PDE remains random, as E $^{σ}$[ u ] is a function of σ , which itself is random. Integrating vertically, we obtain an ini- tial value problem with random parameters for the conditional mean Ekman transport, E $^{σ}$[ T$_{e}$ ]( t ) = ∫ − H 0 E $^{σ}$[ u$_{e}$ ]( z, t ) dz . This conditional mean transport satisfies the following equation:

$$
\frac{d}{dt} \mathbb{E}^{\sigma}[\mathbf{T}_e] = -if \mathbb{E}^{\sigma}[\mathbf{T}_e] + \frac{1}{\rho_o} (\tau_w + \tau_s^0) - if \mathbf{T}_s, \quad (2.16a)
$$

Its solution is given by

$$
\begin{equation}
\begin{split}
\mathbb{E}^{\sigma}[\mathbf{T}_e](t) = e^{-ift}\mathbb{E}^{\sigma}[\mathbf{T}_e](0) &+ \frac{1}{\rho_o}\int_0^t e^{-if(t-r)}(\boldsymbol{\tau}_w + \boldsymbol{\tau}_s^0)(r) dr \\
& - \mathbf{T}_s(1 - e^{-ift}).
\end{split}
\tag{2.16b}
\end{equation}
$$

We recall that the Itô–Stokes drift u$_{s}$ is assumed to be quasistationary (hence T$_{s}$ is also quasistationary), whereas the wind stress τ$_{w}$ and the vertical diffusion coefficient a$_{zz}$ are potentially time dependent and random (and consequently so is τ 0 s = a$_{zz}$∂$_{z}$ u$_{s}$ | z =0 ). Taking the divergence of the previous equation, we derive the following expression for the condi- tional mean Ekman pumping at the lower boundary:

$$
\begin{align}
\mathbb{E}^{\sigma}[w_e](t) ={}& \cos(ft)\mathbb{E}^{\sigma}[w_e](0) + [\cos(ft) + 1]\nabla \cdot \mathbf{T}_s \nonumber \\
& + \frac{1}{\rho_o}\int_0^t \sin[f(t-r)]\nabla \times (\boldsymbol{\tau}_w + \boldsymbol{\tau}_s^0)(r)dr \nonumber \\
& + \cos(ft)\nabla \times [\mathbb{E}^{\sigma}[\mathbf{T}_e](0) + \mathbf{T}_s] \nonumber \\
& + \frac{1}{\rho_o}\int_0^t \sin[f(t-r)]\nabla \cdot (\boldsymbol{\tau}_w + \boldsymbol{\tau}_s^0)(r)dr. \tag{2.17}
\end{align}
$$

The Ekman pumping induced by the unsteady Ekman– Stokes boundary layer depends on both the divergence and curl components of the initial Ekman transport, the steady Itô–Stokes transport, and the time-dependent effective sur- face stress. Note again that compared to the deterministic case, the additional terms are related to the surface Itô–Stokes stress

τ 0 $_{s}$. These formulations of the expectation of Ekman transport and pumping, conditioned on the small-scale correlation tensor, provide an expression for the system’s response in the “average.” They are still random but no longer include Brown- ian terms. These expressions yield solutions for a given realiza- tion of wind stress and waves, highlighting the expected modifications to the Ekman layer due to the spatial inhomoge- neity of the random field, which will be identified in the follow- ing as arising from waves and turbulent transport.

The stationary and nonstationary models described in this section provide a simple stochastic forcing parameterization, based on rigorous derivations, for coarse-resolution models that inadequately resolve the Ekman boundary layer, due to limited vertical resolution. This forcing mechanism facilitates a seamless extension of deterministic codes for ensemble simula- tions, with moderate software adaptation and computational cost. The capability of these models to handle non-Gaussian distributions, extreme events, and the incorporation of different physical processes is discussed in the following sections.

d. Consistent physical parameterization for uncertainty representation

In previous works (Bauer et al. 2020; Resseguier et al. 2021; Brecht et al. 2021; Li et al. 2023a; Tucciarone et al. 2024), noise was parameterized through a spectral decomposition of the (Hilbert–Schmidt) correlation operator σ . Both stationary and time-evolving representations were constructed using techniques such as proper orthogonal decomposition, dy- namic mode decomposition, wavelet basis functions, or self- similarity assumptions, enabling data-driven or model-based noise definitions. From these noise basis functions, the diffu- sion tensor a and the resulting Itô–Stokes drift v$_{s}$ , as defined in Eq. (2.1c), were derived directly. In contrast, this work adopts an inverse approach: Starting from specific physical pa- rameterizations of the vertical diffusion coefficient a$_{zz}$ and the horizontal Stokes drift u$_{s}$ , the goal is to derive the correspond- ing vertical and horizontal components of the unresolved ran- dom flow component σ dB$_{t}$ . Additionally, randomness in the correlation operator σ , acting as an extra random source, is parameterized using a wind-driven process and a steady distri- bution of surface wave directions. For simplicity and without loss of generality, we focus on one-dimensional boundary layer models in the following sections.

1) RANDOM WIND AND WIND STRESS

We adopt here the random wind parameterization pro- posed by McWilliams and Huckle (2006). The total wind u$_{a}$ is decomposed into its time average ⃗ u$_{a}$ and fluctuations u ′ $_{a}$. Each of the fluctuating wind components is modeled as an indepen- dent Ornstein–Uhlenbeck (OU) process:

$$
d\mathbf{u}'_a = -\frac{1}{T_a} \mathbf{u}'_a dt + \sqrt{\frac{2}{T_a}} \Sigma_a d\mathbf{W}_t, \quad (2.18)
$$

where T$_{a}$ denotes the memory time scale and Σ$_{a}$ stands for the standard deviation of u ′ $_{a}$, both of which determine the transient behavior of the OU process. The term W$_{t}$ = W x t + iW y t
combines two independent real-valued standard Brownian motions $W_t^x$ and $W_t^y$, which are both assumed to be independent of the functional cylindrical Brownian motion $B_t$ involved in the unresolved current component. This independence assumption is supported by the separation of atmospheric and oceanic media as well as their differences in time scales. Furthermore, it ensures that the conditional mean Ekman transport solution discussed in the previous section remains valid.

The OU wind fluctuation $u'_a$ (conditioned on the initial value $u'_0$) is a Gaussian process with the following probability distribution:

$$u'_a|u'_0 \sim \mathcal{N}[u'_0 e^{-t/T_a}, \Sigma_a^2(1 - e^{-2t/T_a})]. \quad (2.19)$$

This shows that, over the long term ($t \gg T_a$), the process converges to a centered and stationary distribution with a bounded variance $\Sigma_a^2$, such that $u_a \sim \mathcal{N}(\bar{u}_a, \Sigma_a^2)$.

The ocean surface wind stress $\tau_w$ follows the simplest bulk formula, which is quadratic with respect to the total wind $u_a$:

$$\tau_w = C_D \rho_a |u_a| u_a = \rho_o u_a^2 e^{i\theta_a}, \quad (2.20)$$

where $C_D$ is the air-sea drag coefficient, $\rho_a$ is the air density, $u_a$ is the friction velocity, and $\theta_a$ denotes the angular direction of the wind stress. Let us emphasize that the resulting parameters $\tau_w$, $u_a$, and $\theta_a$ are random processes with bounded stationary variance in the long run, introducing additional randomness into related parameters, as discussed later.

## 2) VERTICAL DIFFUSION COEFFICIENT

For vertical diffusion, we consider the well-established KPP model (Large et al. 1994; McWilliams and Huckle 2006) to represent the vertical eddy viscosity $a_{zz}$ within the turbulent surface boundary layer. This nonlocal parameterization adheres to the predictions of Monin-Obukhov similarity theory, along with extensions discussed in Fox-Kemper et al. (2022). Under the neutral condition, it is formulated as follows:

$$a_{zz}(z) = c_1 u_* h G(\zeta), \quad \zeta = -z/h, \quad h = c_2 \frac{u_*}{f}, \quad (2.21a)$$

where $c_1$ and $c_2$ are constants determining the amplitude and shear of $a_{zz}$, $h$ is the boundary layer depth, and $G$ is a smooth function of the normalized depth $\zeta$. Note that, in the absence of stratification and surface buoyancy flux in this study, the boundary layer depth is constrained by the Ekman layer depth. We adopt a specific version of $G$ proposed by McWilliams and Huckle (2006), which includes an additional regularization term near the ocean surface:

$$G(\zeta) = \zeta(1 - \zeta)^2 + \mathcal{H}(\zeta_0 - \zeta) \frac{(\zeta - \zeta_0)^2}{2\zeta_0}, \quad (2.21b)$$

where $\mathcal{H}$ denotes the Heaviside step function. The regularization term represents an additional mixing process near the surface, possibly accounting for wave breaking and mixing confined to a shallow layer with small thickness $\zeta_0 h$ (McWilliams and Huckle 2006). Since $u_*$ is random, both the diffusion coefficient $a_{zz}$ and

the boundary layer depth $h$ are also random. As previously analyzed, $u_a^2$ has bounded variance, implying that the magnitudes of $\tau_w$ and $a_{zz}$ are similarly bounded. Note that more complex schemes based on the second-moment closure (SMC) (Mellor and Yamada 1982; Umlauf and Burchard 2005; Harcourt 2013) could also be investigated within a similar framework. Unlike KPP, SMC models incorporate memory effects by solving prognostic equations for turbulent kinetic energy and dissipation, allowing for a more gradual adaptation of mixing intensity to external forcing. Exploring their response to random wind forcing in future work could provide insights into turbulence memory effects and differences in vertical mixing compared to KPP.

## 3) STOKES DRIFT AND WAVE STRESS

Consider a steady, monochromatic surface wave with surface elevation, to leading order in wave steepness, expressed as $\eta = \alpha \cos(kx - \omega t)$, where $\alpha$ is the wave amplitude, $k$ is the horizontal wavenumber, and $\omega = (gk)^{1/2}$ is the angular frequency satisfying the deep-water dispersion relation. The horizontal components of the Stokes drift are approximately given by (Phillips 1977)

$$\mathbf{u}_s(z) = U_0 e^{2kz} e^{i\theta_s}, \quad (2.22a)$$

where $U_0 = \omega k \alpha^2$ represents the Stokes drift magnitude and $\theta_s$ is the surface wave propagation direction. This Stokes drift velocity is nonlinear with respect to the wave amplitude and decays exponentially with depth.

In reality, surface gravity waves have a broadband spectrum, resulting in a more complex vertical profile for the Stokes drift (Huang 1971; Jenkins 1989). The KPP model, with a random amplitude due to wind uncertainty, makes the derived wave stress amplitude $\tau_s = a_{zz} \partial_t \mathbf{u}_s$ equally random. To account for uncertainty in the surface wave propagation direction, $\theta_s$ is parameterized using a Gaussian distribution:

$$\theta_s \sim \mathcal{N}(\Theta_s, \Sigma_s^2), \quad (2.22b)$$

where $\Theta_s$ and $\Sigma_s^2$ denote the mean and standard deviation, respectively. Consequently, the wave stress $\tau_s$ becomes random but non-Gaussian due to its nonlinear dependence on $\theta_s$.

## 4) UNRESOLVED RANDOM FLOW

Given the diffusion coefficient $a_{zz}$ and Stokes drift $\mathbf{u}_s$, the vertical and horizontal components of the unresolved noise flow are defined through the following projection formulation:

$$\tilde{\sigma}_z dB_t = \sqrt{2} \sum_n \langle a_{zz}^{1/2}, e_n \rangle e_n d\beta_t^n, \quad (2.23a)$$

$$\tilde{\sigma}_x dB_t = \sqrt{2} \sum_n \langle a_{zz}^{-1/2} U_s, e_n \rangle + i \langle a_{zz}^{-1/2} V_s, e_n \rangle e_n d\beta_t^n, \quad (2.23b)$$

$$\mathbf{U}_s = \int_{-H}^{z} \mathbf{u}_s(\zeta) d\zeta = U_s + iV_s, \quad (2.23c)$$

where $\{e_n\}_n$ is a set of orthogonal basis functions of the real-valued Hilbert space $L^2([0, H], \mathbb{R})$, with inner product $\langle f, g \rangle = \int_{-H}^{0} f(z)g(z)dz$, while $\{\beta_t^n\}_n$ is a set of independent
<figure><img src="image_2.png" /><figcaption>FIG. 1. Illustration of the identifications in Eq. (2.25b), based on the sample variance {$\tilde{a}_{zz} = (1/2)\mathbb{E}[(\sigma_x dB_t)^2]/dt$} and covariance {$\tilde{a}_{xz} = (1/2)\mathbb{E}[(\sigma_x dB_t)(\sigma_z dB_t)]/dt$} of the noise defined in Eq. (2.23), constructed using Dirac delta functions on a Chebyshev grid, as described in appendix B. Note that only the zonal component of the Stokes drift is used in this illustration.</figcaption></figure>

real-valued standard Brownian motions, and $U_s$ denotes the antiderivative of the Stokes drift, reducing to $u_s/(2k)$ for the monochromatic surface wave in Eq. (2.22a). For more complex Stokes drift profiles, $U_s$ can be redefined to enhance accuracy. We emphasize that the independence assumption between the noise terms driving the wind and the ocean currents implies that the Brownian motions {$\beta_t^n$}$_n$ are independent of the pair ($W_t^1, W_t^0$) defined in Eq. (2.18). Note that Eq. (2.23b) is only defined within the support of the function $a_{zz}$, with zero horizontal noise $\sigma_x dB_t$ outside this region. The notation $\tilde{a}$ for the noise variables is used to indicate that the noise is expressed as a projection onto a particular (finite dimensional in practice) basis, with constraints on the associated Itô-Stokes drift and vertical diffusion.

The corresponding diffusion processes, originally defined as $\tilde{a} = (1/2)\mathbb{E}[\tilde{\sigma} dB_t (\tilde{\sigma} dB_t)^T]/dt$, are subsequently expressed as

$$ 
\begin{aligned}
\tilde{a}_{zz} &= \sum_n \langle a_{zz}^{1/2} e_n \rangle^2 e_n^2, \quad \tilde{a}_{xz} \\
&= \sum_n \langle a_{zz}^{1/2}, e_n \rangle (\langle a_{zz}^{-1/2} U_s, e_n \rangle + i \langle a_{zz}^{-1/2} V_s, e_n \rangle) e_n^2.
\end{aligned} 
\quad (2.24) 
$$

Upon using Parseval's theorem, the reconstructed diffusion coefficients are globally consistent with the vertical viscosity and the antiderivative of Stokes drift:

$$ \int_{-H}^{0} \tilde{a}_{zz}(z) dz = \int_{-H}^{0} a_{zz}(z) dz, \quad \int_{-H}^{0} \tilde{a}_{xz}(z) dz = \int_{-H}^{0} U_s(z) dz. \quad (2.25a) $$

If the basis functions {$e_n$} are localized (significant only over small, distinct regions with negligible overlap), Parseval's theorem holds almost pointwise, leading to

$$ \tilde{a}_{zz}(z) \approx a_{zz}(z), \quad \tilde{a}_{xz}(z) \approx U_s(z). \quad (2.25b) $$

Consequently, the resulting Itô-Stokes drift $\tilde{u}_s = \partial_x \tilde{a}_{xz}$ approximates the given Stokes drift $u_s$. This identification is illustrated in Fig. 1, which demonstrates that the sample variance and covariance of the noise approximate the target KPP viscosity and Stokes drift profiles.

While this noise model relies on simplified parameterizations of diffusion (KPP) and Itô-Stokes statistical drift (assumed to match the monochromatic Stokes drift), the distribution of the resulting noise is nontrivial to specify. It exhibits a nonlinear dependence on random parameters, and its sensitivity to small parameter variations is intricate. Furthermore, the noise acts both as an additive source of randomness and as a transport (advection) mechanism, making it inherently non-Gaussian.

### 3. Statistical analyses

This section investigates the statistical properties of the proposed time-dependent stochastic Ekman-Stokes model [Eq. (2.14)]. Large-ensemble simulations are conducted using a Monte Carlo method (Caflisch 1998). The SPDE is solved numerically with a pseudospectral Chebyshev method (Boyd 2001; Trefethen 2000) combined with an implicit time-stepping
<table><caption>TABLE 1. Common parameters used for the simulations.</caption><thead><tr><th>Parameters</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td><i>f</i></td><td>0.73 &times; 10<sup>-4</sup> s<sup>-1</sup></td><td>Coriolis frequency</td></tr><tr><td><i>&rho;</i><sub>o</sub></td><td>1000 kg m<sup>-3</sup></td><td>Water density</td></tr><tr><td><i>&rho;</i><sub>a</sub></td><td>1 kg m<sup>-3</sup></td><td>Air density</td></tr><tr><td><i>C</i><sub>D</sub></td><td>1.3 &times; 10<sup>-3</sup></td><td>Air-sea drag coefficient</td></tr><tr><td><i>u&#773;</i><sub>a</sub></td><td>5 m s<sup>-1</sup></td><td>Mean wind speed</td></tr><tr><td><i>&Sigma;</i><sub>a</sub></td><td>5 m s<sup>-1</sup></td><td>Standard deviation of transient wind</td></tr><tr><td><i>T</i><sub>a</sub></td><td>1 day</td><td>Memory time of fluctuation wind</td></tr><tr><td><i>c</i><sub>1</sub></td><td>0.4</td><td>von Kármán constant</td></tr><tr><td><i>c</i><sub>2</sub></td><td>0.7</td><td>Boundary layer depth parameter</td></tr><tr><td><i>&zeta;</i><sub>0</sub></td><td>0.05</td><td>Normalized depth of shallow layer</td></tr><tr><td><i>&kappa;</i><sub>b</sub></td><td>10<sup>-4</sup> m<sup>2</sup> s<sup>-1</sup></td><td>Background uniform viscosity</td></tr><tr><td><i>&alpha;</i></td><td>0.8 m</td><td>Surface wave amplitude</td></tr><tr><td><i>k</i></td><td>1.05 m<sup>-1</sup></td><td>Wavenumber of Stokes drift</td></tr><tr><td><i>&Theta;</i><sub>s</sub></td><td>0&deg;</td><td>Stokes drift mean direction angle</td></tr><tr><td><i>&Sigma;</i><sub>s</sub></td><td>5&deg;</td><td>Standard deviation of wave angles</td></tr><tr><td><i>H</i></td><td>256 m</td><td>Vertical domain depth</td></tr><tr><td><i>N</i><sub>z</sub></td><td>512</td><td>Number of Chebyshev points</td></tr><tr><td><i>N</i><sub>r</sub></td><td>1000</td><td>Number of random realizations</td></tr><tr><td><i>&Delta;t</i></td><td>30 min</td><td>Time step</td></tr></tbody></table>

scheme (Jentzen and Kloeden 2011). The Markovian wind stress process is integrated using an exact scheme based on its transition distribution [Eq. (2.19)]. Detailed descriptions of these numerical schemes are provided in appendix B. To investigate the influence of stochastic terms, we begin by comparing the ensemble statistics of the proposed model with those of the deterministic Ekman-Stokes model driven by random wind forcing. Next, we conduct parametric studies to assess the effects of random wind and surface wave parameters within this framework. Finally, the model’s ensemble forecasts are assessed using observational data.

### a. Comparison with the deterministic Ekman-Stokes model driven by random wind

The ensemble statistics derived from the proposed parameterized SPDE [Eq. (2.14)], referred to as p-SPDE, are first compared with those of a benchmark model (McWilliams et al. 1997) solely considering the Coriolis-Stokes force, driven by the parameterized PDE, designated as p-PDE:

$$ \partial_t \mathbf{u}_e + \text{if}(\mathbf{u}_e + \mathbf{u}_s) = \partial_z (a_{zz} \partial_z \mathbf{u}_e). \quad (3.1) $$

This comparison is performed under identical boundary conditions as specified in Eq. (2.14b), with both models initialized from a state of rest. The common parameters between both simulations are listed in Table 1. Specifically, we assume the mean directions for both wind velocity and Stokes drift to be zonal. It is noteworthy that the selected surface wave amplitude $\alpha = 0.8$ m and wavelength $\lambda = 60$ m (or wavenumber $k = 2\pi/\lambda \approx 1.05 \text{ m}^{-1}$) proposed by McWilliams et al. (1997) imply a Stokes drift magnitude $U_0 \approx 0.068 \text{ m s}^{-1}$ and a Stokes-layer depth $h_s = 1/(2k) \approx 4.775$ m. Additionally, it is important to mention that the transient wind level in this scenario is set equal to the mean wind level ($\Sigma_a = \bar{u}_a = 5 \text{ m s}^{-1}$) to highlight the rectification of the Ekman layer (McWilliams and Huckle 2006).

Under these parameters, the probability density functions (PDFs) of the surface forcings (wind, wind stress, and the derived wave stress in the p-SPDE) are presented in Fig. 2 after reaching statistically stationary states. These bivariate PDFs are estimated from the ensemble vector fields using a non-parametric Kernel density estimation (KDE) method with Gaussian kernels. As anticipated, Fig. 2a illustrates the Gaussian nature of the wind vector $\mathbf{u}_a$ with uncorrelated components, implying independence in the Gaussian case. However, the nonlinearity of the bulk formula [Eq. (2.20)] results in a non-Gaussian distribution for the wind stress $\tau_w$. Figure 2b reveals that the distribution is skewed toward the left of the mean value. Additionally, Fig. 2c highlights a significant bias in the distribution of the surface wave stress vector $\tau_s^0 = \rho_o a_{zz} \partial_z \mathbf{u}_s|_{z=0}$ derived in the p-SPDE (2.14).

With these parameterized distributions, we proceed to the analysis of the response of the two random models. Figure 3 presents the ensemble mean and spread (mean $\pm$ one standard deviation, indicating variability and uncertainty) of the Ekman velocity components. It also shows the ensemble mean kinetic energy (MKE; representing the energy of the mean solution) and eddy kinetic energy (EKE; representing the variance components) densities, along with the Pearson correlation coefficient (PCC) for the velocity components.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Contour plots of PDF for (a) the wind $u_a$, (b) the wind stress $\tau_w$, and (c) the derived surface wave stress $\tau_s^0$ after 30 days of p-SPDE simulation. In each panel, darker-colored contours represent higher density and + marks indicate the ensemble mean.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Comparison of the (a),(b) ensemble mean (solid lines) and spread (shaded areas) for the Ekman velocity components, (c) the ensemble MKE density (solid lines) and EKE density (dashed lines), as well as (d) the PCC for the Ekman velocity components, using different random models (represented by different colors). Note that these ensemble statistics are averaged over the last 20 days of a 30-day simulation.</figcaption></figure>

Compared to the p-PDE, the p-SPDE considering the additional wave mixing and stochastic transport effects produces a smoother profile of mean Ekman velocities with a higher variability of the realizations (ensemble spread) over depth (Figs. 3a,b). This variability provides an indication on the response of the model with respect to the random terms and gives, hence, a quantification of the model uncertainty. Note that the ensemble spread should be in balance with the model representation error. A too narrow spread does not allow to capture the system’s error while an excessive spread indicates a too high sensitivity. As shown in section 3d, the p-SPDE provides a good error/spread balance.

Additionally, the p-SPDE yields higher MKE and EKE with larger differences between them (Fig. 3c) and stronger correlation between the zonal and meridional components

(Fig. 3d) throughout the depth. These differences are particularly pronounced near the ocean surface. Note that both models exhibit higher EKE than MKE throughout the depth, a negative correlation in the upper 50 m, and a positive correlation in the subsequent 100 m.

Higher-order ensemble statistics of the random models can also be analyzed. The local skewness and kurtosis of the Ekman velocity components over time and depth are shown in Fig. 4. Skewness quantifies the asymmetry of distribution relative to its mean, while (excess) kurtosis characterizes the “tailedness” of the distributions, that is, the prominence of extreme values. Specifically, we compute the sample skewness and kurtosis of a random variable $X$ from an ensemble of $N_r$ realizations $\{X^{(1)}, \dots, X^{(N_r)}\}$. These are denoted by $\hat{s}(X)$ and $\hat{k}(X)$, respectively, and are defined as

<figure><img src="image_5.png" /><figcaption>FIG. 4. Comparison of the ensemble skewness and kurtosis for the Ekman velocity components over time and depth using different random models (grouped by rows). The sample estimators for skewness and kurtosis are defined in Eq. (3.2), with the latter including the standard $-3$ correction. A 1-day low-pass filter is applied to the time series of these statistics at each spatial point.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Evolution comparison after wind stopped for (a),(b) mean Ekman spiral, (c) mean and spread of surface current speed, and (d) mean and spread of Ekman transport magnitude, using different random models. The wind is stopped after 30 days of simulation. The mean Ekman spirals at different days in (a) and (b) are represented by various colors. In each spiral, a darker color indicates a deeper position, with a depth step of 5 m between adjacent points.</figcaption></figure>

$$ \hat{s}(X) = \frac{\hat{m}_3(X)}{\hat{m}_2^{3/2}(X)}, \quad \hat{k}(X) = \frac{\hat{m}_4(X)}{\hat{m}_2^{1/2}(X)} - 3, \\ \hat{m}_p(X) = \frac{1}{N_r} \sum_{i=1}^{N_r} [X^{(i)} - \bar{X}]^p, \qquad (3.2) $$

where $\bar{X} = (1/N_r) \sum_{i=1}^{N_r} X^{(i)}$ is the sample mean and $\hat{m}_p(X)$ denotes the sample pth central moment (with $p > 1$). The "-3" correction in the kurtosis definition (often referred to as excess kurtosis) ensures that the kurtosis of a Gaussian distribution is zero. Positive values indicate distributions with heavier tails (more prone to extreme events) than a normal distribution, while negative values indicate lighter tails. Skewness and kurtosis within the upper 100-m depth are shown in Fig. 4, noting that these normalized statistics become poorly defined in deeper regions due to variance. Compared to the p-PDE, the p-SPDE frequently exhibits greater positive skewness for the zonal component (Figs. 4a,e), greater negative skewness for the meridional component (Figs. 4b,f), and more frequent extremes for both components (Figs. 4c,d,g,h), all predominantly occurring in the upper 50 m.

To further emphasize the differences between the two schemes, the wind is stopped after 30 days of simulation by setting $\tau_w = 0$, while maintaining the last day's coefficient $a_{zz}$ as a steady diffusion together with a stationary Stokes drift $\mathbf{u}_s$. The two random models are then run without surface momentum flux for several more days. The results, Fig. 5, show significant differences in their statistics. The mean Ekman spiral size of the p-PDE (Fig. 5a) quickly shrinks, with the Ekman current profile becoming nearly centered at the origin, whereas the p-SPDE (Fig. 5b) largely retains its mean spiral structure. Only the current speed decreases over time. Figures 5c and 5d show that the mean and spread of the surface current speed, $|\mathbf{u}_s^0|$, as well as the Ekman transport magnitude, $|T_e|$, for both of the two random models rapidly decrease, reaching steady states within 1 week. However, the residual mean values and uncertainties are significantly higher for the p-SPDE than for

the p-PDE. This can be understood from Eqs. (2.10), (2.11) and (2.16b) with $\tau_w = 0$. It emphasizes the role of all the additional terms associated with the Ito-Stokes drift term in maintaining the spiral when the wind stops. In the p-PDE model, the Stokes drift primarily influences the initial and boundary conditions. Its effect diminishes over time and with increasing depth. In contrast, in the p-SPDE model, the Stokes drift enhances vertical diffusion, maintaining a more significant impact throughout the system. These results show the capability of p-SPDE to transition from a wind- and wave-driven model to a wave-driven model when the wind ceases.

b. Wind-driven rectifications

This subsection explores the impact of wind gustiness on the p-SPDE. Sensitivity analyses are performed to assess this effect. For this purpose, various ensembles are simulated over a range of $\Sigma_a$ values while keeping the other parameters in Table 1 invariant. As shown in Figs. 6a and 6b, a greater wind variability reduces the vertical shear of the mean Ekman current, drives deeper circulation, and increases uncertainty. Figure 6c demonstrates that increased wind variability results in larger MKE and EKE throughout the depth, particularly near the surface. Figure 6d shows a weaker negative correlation at the surface but with stronger variations. The correlation support extends with the increase of wind variation, while the correlation peak decreases and shifts deeper in depth. As described in section 2d, increased gustiness results in a higher diffusion coefficient $a_{zz}$ (and hence a higher magnitude of the correlation $\sigma_z$ for the unresolved motions). It also leads to a larger associated surface wave stress $\tau^0 = (a_{zz} \partial_z \mathbf{u}_s)|_{z=0}$, which amplifies both wave mixing and stochastic transport effects in the p-SPDE, as can be seen from Eq. (2.14a) in the vertical diffusion term and the transport random term, which involves both the vertical shear of the Itô-Stokes drift.

The ensemble skewness and kurtosis of the zonal Ekman velocity component with respect to different $\Sigma_a$ values are shown in Fig. 7. More skewed distributions with higher kurtosis are obtained, reflecting an increase in extreme events over
<figure><img src="image_7.png" /><figcaption>FIG. 6. Comparison of the (a),(b) ensemble mean (solid lines) and spread (shaded areas) for the Ekman velocity components, (c) the MKE density (solid lines) and EKE density (dashed lines), and (d) the correlation coefficient for the Ekman velocity components, with respect to different transient wind levels (represented by different colors).</figcaption></figure>

time and near the surface. Similar results were observed for the meridional component (not shown here for concision purpose).

To further quantify the sensitivity across a broader range of $\Sigma_a$ values, we focus on the time-averaged ensemble statistics of diagnostic variables. Figure 8a shows that both global MKE and EKE increase with higher transient wind, with EKE rising more rapidly than MKE, especially when the transient wind exceeds half of the mean wind component. Figure 8b demonstrates that both the mean and variance of the estimated boundary layer depth increase (quasilinearly) with the wind variation level. Figure 8c illustrates that the Ekman transport magnitude $|\mathbf{T}_e|$ grows significantly with wind. Figure 8d indicates that the mean angle of $\mathbf{T}_e$ is approximately 100° (relative to the left of the mean wind direction) with a slight increasing trend, while the uncertainty of this angular distribution rises rapidly with wind.

### c. Surface-wave-induced modifications

In this section, we perform the sensitivity analyses with respect to the wavelength $\lambda$ and the mean angle $\Theta_s$ of the Stokes drift, respectively. The same metrics as in the previous analysis are used here.

First, we examine the impact of wavelength. Figure 9 shows that shorter surface gravity wave results in smoother vertical profiles of the mean Ekman current with higher uncertainty (Figs. 9a,b), higher densities of MKE and EKE (Fig. 9c), and stronger correlation (Fig. 9d), particularly near the surface. Furthermore, as illustrated in Figs. 10a and 10c, both global MKE and EKE, along with the mean and uncertainty of the transport magnitude, decrease with increasing wavelength. Figure 10b shows a slight decreasing trend for the mean and variance of the estimated Ekman layer depth with increasing wavelength. An even subtler decreasing trend is evident in Fig. 10d for the Ekman transport angles. Notably, shorter wavelengths yield a greater magnitude and vertical shear of the Stokes drift [Eq. (2.22a)], thereby amplifying the wave mixing effects in the p-SPDE [Eq. (2.14)].

Now, we examine the sensitivity of the p-SPDE to the mean angular direction of surface waves. Figures 11a and 11c demonstrate that both MKE and EKE, along with the mean and uncertainty of the transport magnitude, reach a maximum when the mean angular direction of the Stokes drift coincides with that of the wind velocity, decreasing symmetrically with higher rotations to either side. Notably, in the aligned case, the derived surface wave stress simply enhances the wind

<figure><img src="image_8.png" /><figcaption>FIG. 7. Comparison of the ensemble kurtosis for the zonal Ekman velocity component over time and depth with respect to different transient wind levels (grouped by columns). The sample estimators for skewness and kurtosis are defined in Eq. (3.2).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 8. (a)–(d) Comparison of the ensemble statistics for diagnostic variables across various transient wind levels. The MKE and EKE are represented by different colors in (a). The mean and uncertainty are depicted using error bars in (b)–(d). These ensemble statistics are averaged over the last 20 days. It is noteworthy that circular statistics (Fisher 1993) are specifically computed in (d).</figcaption></figure>

stress. Unsurprisingly, the response of the Ekman transport angle follows the variation of the surface wave mean directions, as shown in Fig. 11d. The difference between these output and input angles remains nearly invariant.

The p-SPDE model of the Ekman boundary layer has demonstrated its ability to quantify the system’s response to wind stress, turbulent mixing, and wave mixing, particularly in capturing more frequent extreme statistical events. Such a model could be highly valuable for data analysis or as a simple yet physically relevant stochastic forcing parameterization for coarse-resolution models. It would enable, through this, the ensemble simulation with natural variability driven by waves, turbulence, and wind forcing. Its efficiency and accuracy will be further assessed in a future work using references provided by large-eddy simulations. In the following section, we provide an initial assessment of this model with respect to real data.

### d. Comparison with LOTUS observations

This section compares the model outputs with observations from the Long-Term Upper-Ocean Study (LOTUS) experiment (Price et al. 1987), conducted in the western Sargasso Sea (34°N, 70°W) during the summer of 1982. This dataset spans 160 days and includes current profiles measured with

vector-measuring instruments on a stable platform to reduce errors from mooring movement. To isolate the wind-driven component, the geostrophic velocity—assumed constant at 50-m depth—was removed. Daily averages of wind and current data were rotated to align the wind direction with nominal north and then averaged over the 160-day period. However, we should remark that these observations also include uncertainties. Table 1 of Price et al. (1987) provides the mean (time averaged) values at various depths, accompanied by their respective confidence intervals (CIs). Further details about the dataset and its processing are provided in Price and Sundermeyer (1999).

For the model simulations, a mean zonal wind speed of $\bar{u}_a = 6 \text{ m s}^{-1}$, corresponding to the time-averaged wind speed observed in the dataset (Lewis and Belcher 2004), was used as the forcing. The wind stress in the LOTUS observations (Briscoe and Weller 1984) was calculated using bulk aerodynamic formulas based on wind data recorded at 3-m height, adjusted to 10 m using neutral stability assumptions. Under this configuration, the observed zonal and meridional current components represent the downwind and crosswind (with opposite sign) directions, respectively, as described in Table 1 of Price et al. (1987). To introduce uncertainty, wind variability was added in the meridional direction with a standard

<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison of the (a),(b) ensemble mean (solid lines) and spread (shaded areas delimited by dotted colored lines) for the Ekman velocity components, (c) the MKE density (solid lines) and EKE density (dashed lines), and (d) the correlation coefficient for the Ekman velocity components, with respect to different wavelengths.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Comparison of the ensemble statistics for diagnostic variables across various wavelengths.</figcaption></figure>

deviation of $\Sigma(v'_n) = 1 \text{ m s}^{-1}$. The Coriolis parameter at the observation site is $0.836 \times 10^{-4} \text{ s}^{-1}$.

To handle the shallow surface layer, we used a KPP scheme (McWilliams et al. 2012) without the additional regularization term in Eq. (2.21b). This choice was motivated by the sensitivity of the additional mixing term to the near-surface shallow-layer depth parameter $\zeta_0$, which poses challenges for comparison with LOTUS observations starting from 5 m below the surface. To avoid numerical issues, the computational domain was set to start at 1-m depth and extend to 100 m, discretized with 300 Chebyshev points. Other model parameters were consistent with those in Table 1.

Within a stochastic framework, model accuracy is often assessed by verifying ensemble forecasts. Here, we use a simple ensemble reliability criterion: the mean-square error (MSE) of the ensemble-mean forecast should approximately match the mean ensemble variance (VAR; also known as mean-square dispersion), up to a scaling factor dependent on the ensemble size (Stephenson and Doblas-Reyes 2000; Weigel 2012; Fortin et al. 2014). This is expressed as

$$ \text{MSE} := \frac{1}{N_p} \sum_{n=1}^{N_p} (\bar{q}_n - q_n^o)^2 \approx \frac{N_r + 1}{N_r} \left[ \frac{1}{N_p} \sum_{n=1}^{N_p} \hat{\sigma}^2(q_n) \right] := \text{VAR}, \quad (3.3) $$

where $\bar{q}_n = (1/N_r) \sum_{m=1}^{N_r} q_m^{(m)}$ and $\hat{\sigma}^2(q_n) = [1/(N_r - 1)] \sum_{m=1}^{N_r} [q_m^{(m)} - \bar{q}_n]^2$ represent the empirical estimators of the ensemble mean and variance, $q_m^{(m)}$ is the forecast from the $m$th ensemble member at the $n$th point, and $q_n^o$ is the observation. For our application, the observed and forecasted quantities are the Ekman velocities at discrete depths, i.e., $q_n^{(m)} = \mathbf{u}_e^{(m)}(z_n)$ and $q_n^o = \mathbf{u}_e^o(z_n)$. Specifically, we focus on the mean values of the LOTUS observations in this context. To

quantify forecast performance, we use the absolute difference between the MSE and VAR: $\epsilon = |\text{MSE} - \text{VAR}|$. This metric indicates how well the forecast system balances model error representation and ensemble spread dispersion.

We compared this metric for three models: the p-PDE, the p-SPDE, and a p-SPDE without wave mixing terms (denoted as p-SPDE-noWM), described by

$$ d\mathbf{u}_e = [-\text{if}(\mathbf{u}_e + \mathbf{u}_s) + \partial_z(a_{zz}\partial_z\mathbf{u}_e)]dt - (\text{if}\sigma_x + \sigma_z\partial_z\mathbf{u}_e)dB_t. \quad (3.4) $$

This comparison highlights the contributions of stochastic transport and wave mixing in the proposed formulation.

Figures 12a and 12b show the ensemble mean and spread of the different ensemble models, along with the observations at various points. It can be observed that the p-SPDE model, which includes wave mixing effects, produces a broader ensemble spread that encompasses most of the mean observations and overlaps their CIs, compared to the other models. Additionally, Fig. 12c shows that the p-SPDE achieves the lowest difference between MSE and VAR, indicating a better balance between model error representation and ensemble dispersion. After 5–6 days, an almost perfect balance is achieved. The p-SPDE-noWM also improves this metric compared to the p-PDE, highlighting the contribution of stochastic transport terms. The balance reaches equilibrium at a slightly higher level than the p-SPDE. However, it still improves the p-PDE solution.

Finally, similar tests were conducted with reduced surface wave effects (i.e., with a reduced Stokes drift). The results, shown in Figs. 13a and 13b, indicate that the ensemble spread for all models is reduced compared to that observed in the previous case (Fig. 12). In particular, no model is able to encompass the mean observation of $v_e$ or overlap its CI at $z = -5 \text{ m}$. The modification of the p-SPDE near the surface

<figure><img src="image_12.png" /><figcaption>FIG. 11. Comparison of the ensemble statistics for diagnostic variables across various mean angles.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12. Comparison of ensemble forecasts from different random models with LOTUS observations under wave parameters $U_s = 0.068 \ \text{m s}^{-1}$ and $h_s = 4.775 \ \text{m}$ (McWilliams et al. 1997). The ensemble mean and spread are shown by solid lines and shaded areas, respectively. Black error bars indicate the CIs of the observed means (central points), based on the values reported in Table 1 of Price et al. (1987). The forecast reliability is assessed by the normalized absolute difference between the MSE of the ensemble mean relative to the observed mean and the mean intraensemble variance, scaled by the ensemble size. Lower values indicate better forecast performance.</figcaption></figure>

due to wave mixing effects is less pronounced than in the previous case. Figure 13c shows that the ensemble forecast score values of all models increase when surface waves become smaller. Nevertheless, the rank of model performance remains unchanged from the previous case, with the p-SPDE outperforming the other two models. This result further validates the proposed framework.

## 4. Conclusions

In this study, we applied the stochastic framework of location uncertainty to model the upper-ocean Ekman boundary layer and its variability. This formulation couples the contributions of

random fluctuations with the boundary layer dynamics via an uncertainty representation of unresolved motions, fully consistent with established physical parameterizations.

Through numerical investigations, the statistical responses of the proposed time-dependent stochastic Ekman–Stokes model are analyzed. Comparisons are performed with a benchmark model solely driven by the Coriolis–Stokes force. Incorporating wave mixing and stochastic transport compound effects, the stochastic model produces mean Ekman velocity profiles with a stronger vertical gradient for most depths, along with increased uncertainty and higher mean kinetic energy (MKE) and eddy kinetic energy (EKE). Additionally, numerical results display stronger correlations between zonal and meridional

<figure><img src="image_14.png" /><figcaption>FIG. 13. Comparison of ensemble forecasts from different random models with observations under wave parameters $U_s = 0.035 \ \text{m s}^{-1}$ and $h_s = 4 \ \text{m}$ (representing shorter wave effects).</figcaption></figure>
components, with higher-order statistics indicating more skewed distributions and extreme value occurrences, particularly near the surface. The role of the Stokes drift is shown to be strengthened by an enhanced diffusion, maintaining the Ekman spiral when transitioning from a wind- and wave-driven Ekman model to a wave-driven model when the wind ceases.

Sensitivity analyses revealed that increased wind variability reduces vertical shear, deepens circulation, and increases uncertainty and higher-order moments. Similarly, shorter surface waves lead to higher MKE, EKE, and stronger correlations. Both MKE, EKE, and Ekman transport reach their maximum when the mean angular direction of the surface waves coincides with that of the wind velocity. As the angular difference between the wind velocity and surface wave direction increases, these quantities decrease, resulting in a redistribution of current velocity statistics. In terms of sensitivity, this study suggests (in the range of the parameters explored) that the primary source of uncertainty is associated with the wind, followed by the surface waves, with their direction being the final source of uncertainty. The results suggest that the proposed stochastic model could enhance the generalized Ekman model, as previously outlined in the literature, by achieving a better balance between model error representation and ensemble dispersion.

The proposed stochastic model can be applied within the Ekman boundary layer to parameterize random forcing associated with wind stress, turbulence, and wave mixing effects. This approach provides a cost-effective, physically consistent parameterization of unresolved small-scale processes for low-resolution ensemble simulations, which can be utilized for data assimilation purposes or to simulate plausible sets of realizations.

Looking ahead, future research avenues include extending the stochastic Ekman-Stokes model to incorporate stratification effects (Price and Sundermeyer 1999; McWilliams et al. 2009; Gula et al. 2014; McWilliams et al. 2015) and exploring the distinct impacts of wave mixing and stochastic transport on thermal front evolution (Crowe and Taylor 2018, 2019). It is important to note that this study does not account for the feedback effects of surface waves, potentially modulated by upper-ocean random currents, on wind stress, particularly regarding the atmospheric wave-induced turbulent components (Ayet and Chapron 2022). Within the proposed stochastic framework, developing more accurate noise terms to better represent fully coupled ocean-atmosphere Ekman models (Lewis and Belcher 2004) would be a valuable direction for future research. This stochastic formulation can be applied to both the atmosphere and ocean boundary layers, coupled with the parameterized bulk formulation at the air-sea interface, with dependence on the resulting random wind and current.

Additionally, further investigation of the nonlinear stochastic Craik-Leibovich equations through numerical studies, particularly employing large-eddy simulations (LES) to include ocean Langmuir circulation (McWilliams et al. 1997; Harcourt and D'Asaro 2008; McWilliams et al. 2012; Sullivan and McWilliams 2019), could provide deeper insights into the coupled dynamics of wave-current interactions and improve the representation of turbulent processes in ocean models. The

present stochastic modeling approach can serve as an alternative SGS model that accounts for unresolved fluctuations within the LES framework. This approach offers promising prospects for advancing our understanding of turbulent ocean processes and improving predictive capabilities.

*Acknowledgments.* The authors sincerely thank the two anonymous reviewers for their insightful comments and suggestions, which have greatly contributed to improving the manuscript. We also acknowledge the support of the ERC EU Project 856408-STUOD, which made this research possible.

*Data availability statement.* The code to reproduce the simulation data, diagnostics, and figures can be accessed at https://github.com/matlong/SGESM.

## APPENDIX A

### Derivation of Nondimensional Stochastic Momentum Equations

The procedure for deriving the stochastic Ekman-Stokes models begins with scaling the generalized stochastic Craik-Leibovich momentum equations [Eq. (2.1)], assuming constant buoyancy ($b = 0$). To streamline the derivation, we first apply the vector calculus identity, $\nabla(\mathbf{A} \cdot \mathbf{B}) = (\mathbf{A} \cdot \nabla)\mathbf{B} + (\mathbf{B} \cdot \nabla)\mathbf{A} + \mathbf{A} \times (\nabla \times \mathbf{B}) + \mathbf{B} \times (\nabla \times \mathbf{A})$, to Eq. (2.1a). This results in the following form:

$$ 
\begin{aligned}
d\mathbf{v} + (\mathbf{v}dt + \boldsymbol{\sigma}dB_t) \cdot \nabla(\mathbf{v} + \mathbf{v}_s) + f\hat{\mathbf{z}} \times [(\mathbf{v} + \mathbf{v}_s)dt + \boldsymbol{\sigma}dB_t] \\
= -\nabla(pdt + dP_t) + \nabla \cdot [\mathbf{a}\nabla(\mathbf{v} + \mathbf{v}_s)]dt,
\end{aligned} 
\quad (A.1) $$

where $\mathbf{v} = (u, v, w)^T$ and $\nabla = (\partial_x, \partial_y, \partial_z)^T$. The horizontal variables and Coriolis frequency are first scaled as $\mathbf{x} = (x, y)^T = L\hat{\mathbf{x}}$, $\mathbf{u} = (u, v)^T = U\hat{\mathbf{u}}$, and $f = f_0\hat{f}$. Here, the capital letters represent variable scales and hatted variables denote non-dimensional variables. Following Pedlosky (1990), we adopt a change of vertical coordinates to account for stretching in the boundary layer so that $z = H\hat{z}$, $\hat{\zeta} = \hat{z}/\delta_e$, and $\delta_e = H\delta_e$. Here, $\hat{z} \neq O(1)$, $\hat{\zeta} = O(1)$, and $\delta_e$ denotes the Ekman boundary layer thickness. The scales of vertical velocity $w = (UH\delta_e/L)\hat{w}$ and dynamic pressure $p = f_0UL\hat{p}$ are suggested, respectively, by the continuity equation [Eq. (2.1b)] and the geostrophic balance.

The correlation operator is decomposed into horizontal and vertical components, $\boldsymbol{\sigma} = (\boldsymbol{\sigma}_x, \boldsymbol{\sigma}_z)^T$ with $\boldsymbol{\sigma}_x = (\sigma_x, \sigma_y)^T$. The derived diffusion tensor can thus be decomposed into $\mathbf{a}_{xx} = (1/2)\boldsymbol{\sigma}_x\boldsymbol{\sigma}_x^T$, $a_{zz} = (1/2)\sigma_z^2$, and $\mathbf{a}_{xz} = (1/2)\boldsymbol{\sigma}_x\boldsymbol{\sigma}_z$. We propose scaling $\mathbf{a}_{xx}$ and $a_{zz}$ as horizontal and vertical eddy viscosity coefficients, as

$$ \mathbf{a}_{xx} = A_h \hat{\mathbf{a}}_{xx}, \quad a_{zz} = A_v \hat{a}_{zz}, \quad (A.1a) $$

respectively. The scaling for the correlation operator and the cross component of the diffusion tensor follows:

$$ \boldsymbol{\sigma}_x = \sqrt{2}A_h^{1/2}\hat{\boldsymbol{\sigma}}_x, \quad \boldsymbol{\sigma}_z = \sqrt{2}A_v^{1/2}\hat{\boldsymbol{\sigma}}_z, \quad \mathbf{a}_{xz} = A_h^{1/2}A_v^{1/2}\hat{\mathbf{a}}_{xz}. \quad (A.1b) $$
We next recall the Rossby number Ro, horizontal Ekman number $E_h$, vertical Ekman number $E_v$, and the aspect ratio $\gamma$, which are defined by

$$Ro = \frac{U}{f_0 L}, \quad E_h = \frac{A_h}{f_0 L^2}, \quad E_v = \frac{A_v}{f_0 H^2}, \quad \gamma = \frac{H}{L}. \quad (A.2)$$

Using these dimensionless numbers, the Itô–Stokes drift $v_s$ normalized by the characteristic scales of $v$ scales as

$$\tilde{\mathbf{u}}_s := \frac{\mathbf{u}_s}{U} = \frac{E_h}{Ro} \hat{\nabla} \cdot \overline{\mathbf{a}}_{xx} + \frac{E_h^{1/2} E_v^{1/2}}{Ro} \frac{\partial_z \overline{\mathbf{u}}_{xz}}{\delta_e}, \quad (A.3a)$$

$$\tilde{w}_s := \frac{w_s}{W} = \frac{E_h^{1/2}}{Ro} E_v^{1/2} \hat{\nabla} \cdot \overline{\mathbf{a}}_{xz} + \frac{1}{Ro} \frac{E_v}{\delta_e} \frac{\partial_z \overline{a}_{zz}}{\delta_e}. \quad (A.3b)$$

Note that $\hat{\nabla} = (\partial_z \partial_\hat{\gamma})^T$ represents the nondimensional horizontal gradient from here on.

The advection time scale $t = (L/U)\hat{t}$ is imposed as the characteristic time scale. Note that the time scale $\hat{t}$ is not necessarily of $O(1)$. The variance of a (cylindrical) Brownian motion also takes the scale of time, with $B_t = (L/U)^{1/2} \hat{B}_t$. Moreover, the scale of unresolved pressure noise are suggested by the geostrophic balance $f\hat{z} \times \sigma dB_t \approx -\nabla dP_t$, hence $P_t = \sqrt{2}(E_h/Ro)^{1/2} f_0 L^2 \hat{P}_t$.

Substituting the above scalings in Eq. (A.1), the following nondimensional stochastic momentum equations are obtained:

$$\begin{aligned} \mathrm{Ro}[d\hat{\mathbf{u}} + (\hat{\mathbf{u}} \cdot \hat{\nabla} + \hat{w}\partial_z)(\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s)d\hat{t}] \\ + \sqrt{2}\mathrm{Ro}^{1/2} \left( E_h^{1/2} \overline{\sigma}_x \cdot \hat{\nabla} + \frac{E_v^{1/2}}{\delta_e} \overline{\sigma}_z \partial_z \right) (\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s) d\hat{B}_t \\ = -[\hat{f}(\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s)]^\perp + \hat{\nabla} \hat{p}] d\hat{t} - \sqrt{2} \frac{E_h^{1/2}}{\mathrm{Ro}^{1/2}} (\hat{f} \overline{\sigma}_x^\perp d\hat{B}_t + \hat{\nabla} d\hat{P}_t) \\ + E_h \hat{\nabla} \cdot [\overline{\mathbf{a}}_{xx} \overline{\nabla} (\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s)] d\hat{t} + \frac{E_v^{1/2}}{\delta_e} \partial_z \overline{[a}_{zz} \partial_z (\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s)] d\hat{t} \\ + E_h^{1/2} \frac{E_v^{1/2}}{\delta_e} \{\hat{\nabla} \cdot [\overline{\mathbf{a}}_{xz} \partial_z (\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s)] + \partial_z \overline{[\mathbf{a}}_{xz} \overline{\nabla} (\hat{\mathbf{u}} + \tilde{\mathbf{u}}_s)]\} d\hat{t}, \end{aligned} \quad (A.4a)$$

$$\begin{aligned} \gamma^2 \mathrm{Ro}[d\hat{w} + (\hat{\mathbf{u}} \cdot \hat{\nabla} + \hat{w}\partial_z)(\hat{w} + \tilde{w}_s)d\hat{t}] \\ + \sqrt{2}\gamma^2 \mathrm{Ro}^{1/2} \left( E_h^{1/2} \overline{\sigma}_x \cdot \hat{\nabla} + \frac{E_v^{1/2}}{\delta_e} \overline{\sigma}_z \partial_z \right) (\hat{w} + \tilde{w}_s) d\hat{B}_t \\ = -\frac{1}{\delta_e^2} \left( \partial_z \hat{p} d\hat{t} + \sqrt{2} \frac{E_h^{1/2}}{\mathrm{Ro}^{1/2}} \partial_z d\hat{P}_t \right) \\ + \gamma^2 \left( E_h \overline{\nabla} \cdot [\overline{\mathbf{a}}_{xx} \overline{\nabla} (\hat{w} + \tilde{w}_s)] + \frac{E_v}{\delta_e^2} \partial_z \overline{[\hat{\mathbf{a}}}_{zz} \partial_z (\hat{w} + \tilde{w}_s)] \right) \\ + E_h^{1/2} \frac{E_v^{1/2}}{\delta_e} \{ \overline{\nabla} \cdot [\overline{\mathbf{a}}_{xz} \overline{\nabla} (\hat{w} + \tilde{w}_s)] + \partial_z \overline{[\mathbf{a}}_{xz} \overline{\nabla} (\hat{w} + \tilde{w}_s)] \} d\hat{t}, \end{aligned} \quad (A.4b)$$

where $\hat{\mathbf{u}}^\perp = (-\hat{v}, \hat{u})^T$. To further derive the equations of motion in the Ekman boundary layer, we make the following

assumptions regarding the dimensionless numbers (Pedlosky 1990):

$$\mathrm{Ro} \ll 1, \gamma \ll 1, \delta_e \sim E_v^{1/2}, E_h^{1/2} \sim \mathrm{Ro}. \quad (A.5)$$

Using these assumptions, the previous nondimensional momentum equations reduce to

$$\mathrm{Ro} d\hat{\mathbf{u}} = \left\{ -\hat{f}(\hat{\mathbf{u}} + \partial_z \overline{\mathbf{a}}_{xz})^\perp - \hat{\nabla} \hat{p} + \partial_z \overline{[a}_{zz} \partial_z (\hat{\mathbf{u}} + \partial_z \overline{\mathbf{a}}_{xz})} \right. \\ \left. + \mathcal{O}(\mathrm{Ro})\right\} d\hat{t} - \sqrt{2} \mathrm{Ro}^{1/2} \{[\sigma_z \partial_z (\hat{\mathbf{u}} + \partial_z \overline{\mathbf{a}}_{xz}) + \hat{f} \overline{\sigma}_x^\perp] d\hat{B}_t + \hat{\nabla} d\hat{P}_t\} + \mathcal{O}(\mathrm{Ro}) d\hat{B}_t, \quad (A.6a)$$

$$\gamma^2 \mathrm{Ro} d\hat{w} = \left[ -\frac{1}{\delta_e^2} \partial_z \hat{p} + \mathcal{O}(\gamma^2) \right] d\hat{t} - \sqrt{2} \frac{\mathrm{Ro}^{1/2}}{\delta_e^2} \partial_z d\hat{P}_t \\ + \mathcal{O}(\gamma^2) d\hat{B}_t. \quad (A.6b)$$

Subsequently, a steady Ekman–Stokes model and a time-dependent random model are derived in section 2, based on distinct nondimensional time variables $\hat{t}$.

## APPENDIX B

### Numerical Schemes

A pseudospectral Chebyshev method (Boyd 2001) is employed for the vertical discretization of the boundary layer models. The so-called Chebyshev points, denoted as $\hat{z}_j = \cos(j\pi/N_z)$, $j = 0, ..., N_z$, are depicted in Fig. B1. Geometrically, these points represent the projections on [-1, 1] of equispaced points on the upper half of the unit circle. Hence, the projection nodes are denser near the two boundaries than in the midregions, which proves advantageous for boundary layer problems and for representing the Stokes drift near the ocean surface. The grid points of the vertical domain depth $H$ are constructed using a linear transformation, $z_j = -H(\hat{z}_j + 1)/2$. Then, the Chebyshev spectral derivative and integration operators are built on these points to solve a (S)PDE for nonperiodic functions with high accuracy.

<figure><img src="image_15.png" /><figcaption>FIG. B1. Geometrical visualization of Chebyshev points (red points) as the projections on [-1, 1] of equally spaced points on the upper half of the unit circle.</figcaption></figure>
A linear-implicit Euler scheme (Jentzen and Kloeden 2011, chapter 8) is used for the time stepping of the SPDE (2.14), namely,

$$ \mathbf{u}_j^{n+1} = (I_{j,k} - A_{j,k} \Delta t)^{-1} [\mathbf{u}_k^n + F_k^n(\mathbf{u}_s) \Delta t + G_k^n(\mathbf{u}^n, \mathbf{u}_s) \Delta B]. \quad (B.1) $$

Here, $\mathbf{u}_j^n$ denotes the discrete value of the Ekman current velocity $\mathbf{u}_e$ on the node $z_j$ at time $t_n$ and $\mathbf{A} = -\mathbf{f}\mathbf{I} + \mathbf{D}\mathbf{a}_{zz}\mathbf{D}$ represents the discrete linear operator of the SPDE, where $\mathbf{I}$ denotes the identity matrix and $\mathbf{D}$ represents the Chebyshev differentiation matrix (Trefethen 2000, chapter 6). The terms $F$ and $G$ correspond to the right-hand side forcing term of Eq. (2.14a).

Additionally, we consider an exact scheme for the wind using the Markovian property of the Ornstein-Uhlenbeck process [Eq. (2.18)], expressed as

$$ \mathbf{u}_a^{n+1} = \bar{\mathbf{u}}_a + (\mathbf{u}_a^n - \bar{\mathbf{u}}_a)e^{-\Delta t/T_a} + \sum_a \sqrt{1 - e^{-2\Delta t/T_a}} \xi^n. \quad (B.2) $$

Here, $\xi^n = \xi_x^n + i\xi_y^n$ where $\xi_x^n$ and $\xi_y^n$ are independent random variables following a standard normal distribution.

## REFERENCES

Ayet, A., and B. Chapron, 2022: The dynamical coupling of wind-waves and atmospheric turbulence: A review of theoretical and phenomenological models. *Bound.-Layer Meteor.*, **183**, 1–33, https://doi.org/10.1007/s10546-021-0066-6.

Bauer, W., P. Chandramouli, B. Chapron, L. Li, and E. Mémin, 2020: Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: A stochastic approach. *J. Phys. Oceanogr.*, **50**, 983–1003, https://doi.org/10.1175/JPO-D-19-0164.1.

Bender, C. M., and S. A. Orszag, 1978: *Advanced Mathematical Methods for Scientists and Engineers*. McGraw-Hill, 593 pp.

Boyd, J. P., 2001: *Chebyshev and Fourier Spectral Methods*. 2nd ed. Dover Publications, 668 pp.

Brecht, R., L. Li, W. Bauer, and E. Mémin, 2021: Rotating shallow water flow under location uncertainty with a structure-preserving discretization. *J. Adv. Model. Earth Syst.*, **13**, e2021MS002492, https://doi.org/10.1029/2021MS002492.

Briscoe, M. G., and R. A. Weller, 1984: Preliminary results from the long-term upper-ocean study (LOTUS). *Dyn. Atmos. Oceans*, **8**, 243–265, https://doi.org/10.1016/0377-0265(84)90011-3.

Caflisch, R. E., 1998: Monte Carlo and quasi-Monte Carlo methods. *Acta Numer.*, **7**, 1–49, https://doi.org/10.1017/S0962492900002804.

Chapron, B., P. Dérrian, E. Mémin, and V. Resseguier, 2018: Large-scale flows under location uncertainty: A consistent stochastic framework. *Quart. J. Roy. Meteor. Soc.*, **144**, 251–260, https://doi.org/10.1002/qj.3198.

Chu, P. C., 2008: Probability distribution function of the upper equatorial Pacific current speeds. *Geophys. Res. Lett.*, **35**, L12606, https://doi.org/10.1029/2008GL033669.

Craik, A. D. D., and S. Leibovich, 1976: A rational model for Langmuir circulations. *J. Fluid Mech.*, **73**, 401–426, https://doi.org/10.1017/S0022112076001420.

Crowe, M. N., and J. R. Taylor, 2018: The evolution of a front in turbulent thermal wind balance. Part 1. Theory. *J. Fluid Mech.*, **850**, 179–211, https://doi.org/10.1017/jfm.2018.448.

—, and —, 2019: The evolution of a front in turbulent thermal wind balance. Part 2. Numerical simulations. *J. Fluid Mech.*, **880**, 326–352, https://doi.org/10.1017/jfm.2019.688.

Da Prato, G., and J. Zabczyk, 2014: *Stochastic Equations in Infinite Dimensions*. 2nd ed. Cambridge University Press, 493 pp.

Debusssche, A., B. Hug, and E. Mémin, 2023: A consistent stochastic large-scale representation of the Navier-Stokes equations. *J. Math. Fluid Mech.*, **25**, 19, https://doi.org/10.1007/s00021-023-00764-0.

Ekman, V. W., 1905: On the influence of the Earth's rotation on ocean-currents. *Ark. Mat. Astron. Fys.*, **2**, 1–52.

Fisher, N. I., 1993: *Statistical Analysis of Circular Data*. Cambridge University Press, 277 pp.

Fortin, V., M. Abaza, F. Ancil, and R. Turcotte, 2014: Why should ensemble spread match the RMSE of the ensemble mean? *J. Hydrometeor.*, **15**, 1708–1713, https://doi.org/10.1175/JHM-D-14-0008.1.

Fox-Kemper, B., L. Johnson, and F. Qiao, 2022: Ocean near-surface layers. *Ocean Mixing*, M. Meredith and A. N. Garabato, Eds., Elsevier, 65–94, https://doi.org/10.1016/B978-0-12-821512-8.00011-6.

Gula, J., M. J. Molemaker, and J. C. McWilliams, 2014: Submeso-scale cold filaments in the Gulf Stream. *J. Phys. Oceanogr.*, **44**, 2617–2643, https://doi.org/10.1175/JPO-D-14-0029.1.

Harcourt, R. R., 2013: A second-moment closure model of Langmuir turbulence. *J. Phys. Oceanogr.*, **43**, 673–697, https://doi.org/10.1175/JPO-D-12-0105.1.

—, and E. A. D'Asaro, 2008: Large-eddy simulation of Langmuir turbulence in pure wind seas. *J. Phys. Oceanogr.*, **38**, 1542–1562, https://doi.org/10.1175/2007JPO3842.1.

Higgins, C., J. Vanneste, and T. S. van den Bremer, 2020: Unsteady Ekman-Stokes dynamics: Implications for surface wave-induced drift of floating marine litter. *Geophys. Res. Lett.*, **47**, e2020GL089189, https://doi.org/10.1029/2020GL089189.

Huang, N. E., 1971: Derivation of Stokes drift for a deep-water random gravity wave field. *Deep-Sea. Res.*, **18**, 255–259, https://doi.org/10.1016/0011-7471(71)90115-X.

—, 1979: On surface drift currents in the ocean. *J. Fluid Mech.*, **91**, 191–208, https://doi.org/10.1017/S0022112079000112.

Jenkins, A. D., 1986: A theory for steady and variable wind-and-wave-induced currents. *J. Phys. Oceanogr.*, **16**, 1370–1377, https://doi.org/10.1175/1520-0485(1986)016<1370:ATFSAV>2.0.CO;2.

—, 1989: The use of a wave prediction model for driving a near-surface current model. *Dtsch. Hydrogr. Z.*, **42**, 133–149, https://doi.org/10.1007/BF02226291.

Jentzen, A., and P. E. Kloeden, 2011: *Taylor Approximations for Stochastic Partial Differential Equations*. Society for Industrial and Applied Mathematics, 12 pp., https://doi.org/10.1137/1.9781611972016.

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

Le Gall, J.-F., 2016: *Brownian Motion, Martingales, and Stochastic Calculus*. Graduate Texts in Mathematics, Vol. 274, Springer International Publishing, 273 pp.

Leibovich, S., 1980: On wave-current interaction theories of Langmuir circulations. *J. Fluid Mech.*, **99**, 715–724, https://doi.org/10.1017/S0022112080000857.

Lewis, D. M., and S. E. Belcher, 2004: Time-dependent, coupled, Ekman boundary layer solutions incorporating Stokes drift.
Dyn. Atmos. Oceans, **37**, 313–351, https://doi.org/10.1016/j.
dynatmocoe.2003.11.001.

Li, L., B. Deremble, N. Lahaye, and E. Mémin, 2023a: Stochastic data-driven parameterization of unresolved eddy effects in a baroclinic quasi-geostrophic model. *J. Adv. Model. Earth Syst.*, **15**, e2022MS003297, https://doi.org/10.1029/2022MS003297.

—, E. Mémin, and G. Tissot, 2023b: Stochastic parameterization with dynamic mode decomposition. *Stochastic Transport in Upper Ocean Dynamics*, B. Chapron et al., Eds., Mathematics of Planet Earth, Vol. 10, Springer, 179–193, https://doi.org/10. 1007/978-3-031-18988-3_11.

McWilliams, J. C., and P. P. Sullivan, 2000: Vertical mixing by Langmuir circulations. *Spill Sci. Technol. Bull.*, **6**, 225–237, https://doi.org/10.1016/S1353-2561(01)00041-X.

—, and E. Huckle, 2006: Ekman layer rectification. *J. Phys. Oceanogr.*, **36**, 1646–1659, https://doi.org/10.1175/JPO2912.1.

—, P. P. Sullivan, and C.-H. Moeng, 1997: Langmuir turbulence in the ocean. *J. Fluid Mech.*, **334**, 1–30, https://doi.org/10. 1017/S0022112096004375.

—, E. Huckle, and A. F. Shchepetkin, 2009: Buoyancy effects in a stratified Ekman layer. *J. Phys. Oceanogr.*, **39**, 2581–2599, https://doi.org/10.1175/2009JPO4130.1.

—, —, J.-H. Liang, and P. P. Sullivan, 2012: The wavy Ekman layer: Langmuir circulations, breaking waves, and Reynolds stress. *J. Phys. Oceanogr.*, **42**, 1793–1816, https://doi.org/10. 1175/JPO-D-12-07.1.

—, J. Gula, M. J. Molemaker, L. Renault, and A. F. Shchepetkin, 2015: Filament frontogenesis by boundary layer turbulence. *J. Phys. Oceanogr.*, **45**, 1988–2005, https://doi.org/10.1175/JPO- D-14-0211.1.

Mellor, G. L., and T. Yamada, 1982: Development of a turbulence closure model for geophysical fluid problems. *Rev. Geophys.*, **20**, 851–875, https://doi.org/10.1029/RG020i004p00851.

Mémin, E., 2014: Fluid flow dynamics under location uncertainty. *Geophys. Astrophys. Fluid Dyn.*, **108**, 119–146, https://doi.org/ 10.1080/03091929.2013.836190.

—, L. Li, N. Lahaye, G. Tissot, and B. Chapron, 2024: Linear wave solutions of a stochastic shallow water model. *Stochastic Transport in Upper Ocean Dynamics II*, B. Chapron et al., Eds., Mathematics of Planet Earth, Vol. 11, Springer, 223– 245, https://doi.org/10.1007/978-3-031-40094-0_10.

Pedlosky, J., 1990: *Geophysical Fluid Dynamics*. 2nd ed. Springer,
710 pp.

Phillips, O. M., 1977: *The Dynamics of the Upper Ocean*. Cambridge
University Press, 336 pp.

Price, J. F., and M. A. Sundermeyer, 1999: Stratified Ekman layers. *J. Geophys. Res.*, **104**, 20467–20494, https://doi.org/10.1029/1999JC900164.

—, R. A. Weller, and R. R. Schudlich, 1987: Wind-driven ocean currents and Ekman transport. *Science*, **238**, 1534–1538, https://doi.org/10.1126/science.238.4833.1534.

Reeks, M. W., 1983: The transport of discrete particles in inhomo-
geneous turbulence. *J. Aerosol Sci.*, **14**, 729–739, https://doi.
org/10.1016/0021-8502(83)90055-1.

Resseguier, V., E. Mémin, and B. Chapron, 2017: Geophysical flows under location uncertainty, Part I, II & III. *Geophys. Astrophys. Fluid Dyn.*, **111**, 149–227.

—, L. Li, G. Jouan, P. Dérian, E. Mémin, and B. Chapron, 2021: New trends in ensemble forecast strategy: Uncertainty quantification for coarse-grid computational fluid dynamics. *Arch. Comput. Methods Eng.*, **28**, 215–261, https://doi.org/10.1007/s11831-020-09437-x.

Stephenson, D. B., and F. J. Doblas-Reyes, 2000: Statistical methods for interpreting Monte Carlo ensemble forecasts. *Tellus*, **52A**, 300–322, https://doi.org/10.3402/tellusa.v52i3.12267.

Sullivan, P. P., and J. C. McWilliams, 2019: Langmuir turbulence and filament frontogenesis in the oceanic surface boundary layer. *J. Fluid Mech.*, **879**, 512–553, https://doi.org/10.1017/jfm.2019.655.

Trefethen, L. N., 2000: *Spectral Methods in MATLAB*. Society for Industrial and Applied Mathematics, 165 pp.

Tucciarone, F. L., E. Mémin, and L. Li, 2024: Data driven stochastic primitive equations with dynamic modes decomposition. *Stochastic Transport in Upper Ocean Dynamics II*, B. Chapron et al., Eds., Mathematics of Planet Earth, Vol. 11, Springer, 321–336, https://doi.org/10.1007/978-3-031-40094-0_15.

Umlauf, L., and H. Burchard, 2005: Second-order turbulence closure models for geophysical boundary layers. A review of recent work. *Cont. Shelf Res.*, **25**, 795–827, https://doi.org/10.1016/j.csr. 2004.08.004.

Vallis, G. K., 2017: *Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. 2nd ed. Cambridge University Press, 946 pp.

Weigel, A. P., 2012: *Ensemble Forecasts*. John Wiley and Sons,
141–166.

Wenegrat, J. O., and M. J. McPhaden, 2016a: A simple analytical model of the diurnal Ekman layer. *J. Phys. Oceanogr.*, **46**, 2877–2894, https://doi.org/10.1175/JPO-D-16-0031.1.

—, and —, 2016b: Wind, waves, and fronts: Frictional effects in a generalized Ekman model. *J. Phys. Oceanogr.*, **46**, 371–394, https://doi.org/10.1175/JPO-D-15-0162.1.

Xu, Z., and A. J. Bowen, 1994: Wave- and wind-driven flow in water of finite depth. *J. Phys. Oceanogr.*, **24**, 1850–1866, https://doi.org/10.1175/1520-0485(1994)024<1850:WAWDFI>2.0.CO;2.