# Stochastic Ekman-Stokes Model with Statistical Analyses

Long Li

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

## To cite this version:

Long Li. Stochastic Ekman-Stokes Model with Statistical Analyses. STUOD 2024 - 5th Stochastic Transport in Upper Ocean Dynamics Annual Workshop, Sep 2024, Rennes, France. hal-04957876

HAL Id: hal-04957876

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

Submitted on 20 Feb 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stochastic Ekman-Stokes Model with Statistical Analyses

## 5th Stochastic Transport in Upper Ocean Dynamics (STUOD) Annual Workshop

Long Li¹, Etienne Mémin¹, Bertrand Chapron¹,²

¹ Odyssey Team, Inria Rennes, France

² LOPS, Ifremer Plouzané, France

September 24, 2024, Inria, Rennes, France
# 1. Introduction

2. Stochastic Formulation

3. Uncertainty Representation

4. Statistical Analyses

5. Conclusion
Classical Ekman surface boundary layer [Ekman, 1905]

For steady, linear flow, uniform density and viscosity, driven by surface wind stress:

$$
\mathbf{f} \mathbf{u}_e = \alpha \partial_{zz}^2 \mathbf{u}_e, \quad \rho_o \alpha \partial_z \mathbf{u}_e |_{z=0} = \tau, \quad \mathbf{u}_e |_{z \to -\infty} = 0
$$

- Explicit solution – Ekman spiral:

$$
\left\{
\begin{aligned}
u_e &= \frac{\sqrt{2}}{\rho_o f h_e} e^{z/\hbar e} \left( \tau^x \cos\left(\frac{z}{h_e} - \frac{\pi}{4}\right) - \tau^y \sin\left(\frac{z}{h_e} - \frac{\pi}{4}\right) \right) \\
v_e &= \frac{\sqrt{2}}{\rho_o f h_e} e^{z/\hbar e} \left( \tau^x \sin\left(\frac{z}{h_e} - \frac{\pi}{4}\right) + \tau^y \cos\left(\frac{z}{h_e} - \frac{\pi}{4}\right) \right)
\end{aligned}
\right.
$$

(Ekman layer depth: $h_e = \sqrt{2\alpha/f}$)

- Ekman transport and Ekman vertical velocity:

$$
\mathbf{T}_e = \int \mathbf{u}_e \, \mathrm{d}z = \frac{1}{\rho_o f} \tau^{\perp}, \quad w_e = \int \nabla \cdot \mathbf{u}_e \, \mathrm{d}z = \frac{1}{\rho_o f} \nabla \times \tau
$$

* A beautiful theory, but generally difficult to observe due to measurement challenges and the influence of phenomena not accounted for, such as nonuniform eddy viscosity, inertial and gravity waves, stratification, and convection.

<figure><img src="image_3.png" /></figure>
* **Motivation:** Develop a stochastic formulation for the upper ocean boundary layer that captures interactions between wind, waves, and currents while accounting for unresolved motion uncertainties.

* **Approach:** Derive the boundary layer model through proper scaling of the generalized stochastic Kraik-Leibovich equations, and constrain the uncertainty representation with established physical parameterizations.

* **Reference:** Li, L., Mémin, E., and Chapron, B. (2024). A generalized stochastic formulation of the Ekman-Stokes model with statistical analyses. Preprint. https://hal.science/hal-04672188.
1. Introduction

2. Stochastic Formulation

3. Uncertainty Representation

4. Statistical Analyses

5. Conclusion
## Stochastic Boussinesq momentum equations

$$d\mathbf{v} + (\mathbf{v}^* \cdot \nabla \mathbf{v} - \nabla \cdot (\mathbf{a}\nabla \mathbf{v}) + \mathbf{f}\hat{\mathbf{z}} \times \mathbf{v}) dt + (\boldsymbol{\sigma} \cdot \nabla \mathbf{v} + \mathbf{f}\hat{\mathbf{z}} \times \boldsymbol{\sigma}) dB_t = (b\hat{\mathbf{z}} - \nabla p) dt - \nabla dP_t$$

$$\boldsymbol{a} = \frac{1}{2} \boldsymbol{\sigma} \boldsymbol{\sigma}^T, \quad \mathbf{v}^* = \mathbf{v} - \mathbf{v}_s, \quad \mathbf{v}_s = \nabla \cdot \boldsymbol{a}$$

$$\nabla \cdot \mathbf{v}^* = 0, \quad \nabla \cdot \boldsymbol{\sigma} = 0$$

$\boldsymbol{\sigma} dB_t$: unresolved flow component, $\boldsymbol{a}$: diffusion coefficient, $\mathbf{v}^*$: effective advection velocity, $\mathbf{v}_s$: Itô-Stokes drift, $dP_t$: unresolved pressure component

Applying a change of variable from $\mathbf{v}$ to $\mathbf{v}^*$ and a classical vector calculus identity, and assuming a quasi-stationary $\mathbf{v}_s \Rightarrow$

## Generalized stochastic Craik-Leibovich momentum equations

$$d\mathbf{v}^* + ((\mathbf{v}^* - \mathbf{v}_s) \cdot \nabla \mathbf{v}^* - \nabla \cdot (\mathbf{a}\nabla \mathbf{v}^*) + \mathbf{f}\hat{\mathbf{z}} \times \mathbf{v}^* - b\hat{\mathbf{z}} + \nabla p) dt + (\boldsymbol{\sigma} \cdot \nabla \mathbf{v}^* + \mathbf{f}\hat{\mathbf{z}} \times \boldsymbol{\sigma}) dB_t + \nabla dP_t = (-\underbrace{\mathbf{f}\hat{\mathbf{z}} \times \mathbf{v}_s}_{\text{Coriolis Stokes}} - \underbrace{\nabla (\mathbf{v}^* \cdot \mathbf{v}_s)}_{\text{pressure correction}} + \underbrace{\mathbf{v}_s \times (\nabla \times \mathbf{v}^*)}_{\text{Craik-Leibovich vortex force}} - \underbrace{(\nabla \times \mathbf{v}_s) \times \mathbf{v}^*}_{\text{Itô-Stokes force}} + \underbrace{\nabla \cdot (\mathbf{a}\nabla \mathbf{v}_s)}_{\text{Itô-Stokes diffusive mixing}}) dt - \underbrace{\boldsymbol{\sigma} \cdot \nabla \mathbf{v}_s dB_t}_{\text{Itô-Stokes advective mixing}}$$
# Scaling for the Ekman boundary layer

In this study, we do not account for stratification effects within the boundary layer ($b = 0$).

* Noise scalings derived from eddy viscosity:

$$ 
\begin{aligned}
\boldsymbol{a}_{xx} &= \frac{1}{2} \sigma_x \sigma_x^T \sim A_h, \quad a_{zz} = \frac{1}{2} \sigma_z^2 \sim A_v \\
\Rightarrow \sigma_x &\sim \sqrt{2} A_h^{1/2}, \quad \sigma_z \sim \sqrt{2} A_v^{1/2} \quad \boldsymbol{a}_{xz} \sim A_h^{1/2} A_v^{1/2}
\end{aligned}
$$

* Assumptions on dimensionless numbers:

$$ 
\begin{gathered}
\mathrm{Ro} = \frac{U}{f_0 L} \ll 1, \quad \mathrm{E}_h = \frac{A_h}{f_0 L^2} \sim \mathrm{Ro}^2, \quad \mathrm{E}_v = \frac{A_v}{f_0 H^2} \sim \delta_e^2, \quad \delta = \frac{H}{L} \ll 1 \\
(E_h / \mathrm{Ro} = 1 / \mathrm{Re} \sim \mathrm{Ro})
\end{gathered}
$$

## Nondimensional momentum equations

$$ 
\begin{aligned}
\mathrm{Ro} \, \mathrm{d} \boldsymbol{u} &= \left( -f(\boldsymbol{u} + \partial_z \boldsymbol{a}_{xz})^{\perp} - \nabla p + \partial_z (\boldsymbol{a}_{zz} \partial_z (\boldsymbol{u} + \partial_z \boldsymbol{a}_{xz})) + \mathcal{O}(\mathrm{Ro}) \right) \mathrm{d}t \\
&\quad - \sqrt{2} \mathrm{Ro}^{1/2} \left( \left( \partial_z \partial_z (\boldsymbol{u} + \partial_z \boldsymbol{a}_{xz}) + f \boldsymbol{\sigma}_x^{\perp} \right) \mathrm{d}B_t + \partial_x \mathrm{d}P_t \right) + \mathcal{O}(\mathrm{Ro}) \mathrm{d}B_t \\
\delta^2 \mathrm{Ro} \, \mathrm{d} \boldsymbol{w} &= \left( -\frac{1}{\delta_e^2} \partial_z p + \mathcal{O}(\delta^2) \right) \mathrm{d}t - \sqrt{2} \frac{\mathrm{Ro}^{1/2}}{\delta_e^2} \partial_z \mathrm{d}P_t + \mathcal{O}(\delta^2) \mathrm{d}B_t
\end{aligned}
$$
Generalized steady Ekman-Stokes model

* An intermediate time scale $t = O(1)$ gives the prognostic equations of O(1):

$$f(\mathbf{u}^* + \mathbf{u}_s)^\perp = -\nabla p + \partial_z(a_{zz}\partial_z(\mathbf{u}^* + \mathbf{u}_s)), \quad \mathbf{u}_s = \partial_z\mathbf{a}_{xz}, \quad \partial_z p = 0$$

* Separating the geostrophic and ageostrophic from $\mathbf{u}^* = \mathbf{u}_g + \mathbf{u}_e$, restoring dimensions, applying boundary conditions ⇒ Ekman-Stokes model in complex notation:

$$\begin{cases} \partial_z(a_{zz}\partial_z\mathbf{u}_e) - if\mathbf{u}_e = if\mathbf{u}_s - \partial_z(a_{zz}\partial_z\mathbf{u}_s) \\ \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 \end{cases}$$

* WKB approximation with variation of parameters for Ekman current stress $\tau_e = \rho_o a_{zz} \partial_z \mathbf{u}_e$ [Wenegrat and McPhaden, 2016]:

$$\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\tau_s - a_{zz}\partial_{zz}^2\tau_s)(\zeta) d\zeta$$

$$\tau_s = \rho_o a_{zz} \partial_z \mathbf{u}_s, \quad \theta(z) = \sqrt{if} \int_{-H}^{z} a_{zz}^{-1/2}(\zeta) d\zeta, \quad G: \text{a symmetric Green's function}$$

* Modified integral properties ($T = \int \mathbf{u} \, dz, w = \int \nabla \cdot \mathbf{u} \, dz$):

$$\mathbf{T}_e = -\frac{i}{\rho_o f}(\tau_w + \tau_s^0) - \mathbf{T}_s, \quad w_e(-H) = \frac{1}{\rho_o f} \nabla \times (\tau_w + \tau_s^0) - \nabla \cdot \mathbf{T}_s$$
Time-dependent stochastic Ekman-Stokes model

* Considering a fast time scale by a change of time coordinates [Crowe and Taylor, 2018], $\tau = t/R_o \Rightarrow B_\tau = B_t/R_o^{1/2}$, separating geostrophic and ageostrophic $\Rightarrow$ SPDE:

$$ \left\{ \begin{array}{l} d\mathbf{u}_e = \left( -\text{if}(\mathbf{u}_e + \mathbf{u}_s) + \partial_z(a_{zz}\partial_z(\mathbf{u}_e + \mathbf{u}_s)) \right) dt - (\text{if}\sigma_x + \sigma_z\partial_z(\mathbf{u}_e + \mathbf{u}_s)) dB_t \\ \qquad \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 \end{array} \right. $$

* Statistical properties of the SPDE were investigated numerically

* Conditional mean transport ($\mathbb{E}^\sigma[\bullet] := \mathbb{E}[\bullet|\sigma]$, $\sigma$ is random and independent to $B_t$):

$$ \begin{gathered} \partial_t \mathbb{E}^\sigma[\mathbf{u}_e] = -\text{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)) \\ \Rightarrow \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)}(\tau_w + \tau_s^0)(r) dr - \mathbf{T}_s(1 - e^{-ift}) \end{gathered} $$
1. Introduction

2. Stochastic Formulation

**3. Uncertainty Representation**

4. Statistical Analyses

5. Conclusion
* K-profile parameterization (KPP) for vertical eddy viscosity [McWilliams and Huckle, 2006] :

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

$$G(\zeta) = \zeta(1-\zeta)^2 + H(\zeta_0 - \zeta) \frac{(\zeta - \zeta_0)^2}{2\zeta_0}$$

$u_*$: friction velocity, $h$: boundary layer depth, $H$: Heaviside step function

* Stokes drift corresponding to a monochromatic deep-water wave [Phillips, 1977] :

$$\eta = \alpha \cos(kx - \omega t), \quad \omega = (gk)^{1/2}$$

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

$\eta$: surface elevation, $\alpha$: wave amplitude, $k$: horizontal wavenumber, $\omega$: angular frequency, $U_0$: Stokes magnitude, $\theta_s$: wave propagation direction

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

<figure><img src="image_5.png" /></figure>
* Markovian wind process [McWilliams and Huckle, 2006] :

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

* Random wind stress (upper boundary) via bulk formula :

$$ \tau_w = C_D \rho_a |\mathbf{u}_a| \mathbf{u}_a = \rho_o u_*^2 e^{i\theta_*} \\ \Rightarrow u_* \text{ random} \Rightarrow a_{zz} \text{ random (via KPP)} $$

* Random surface wave stress (derived in SPDE):

$$ \tau_s^0 = a_{zz} \partial_z \mathbf{u}_s |_{z=0} = a_{zz}(0) U_0 e^{i\theta_s}, \quad \theta_s \sim \mathcal{N}(\Theta_s, \Sigma_s^2) $$

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

<figure><img src="image_7.png" /></figure>
# Noise representation

**Objective:** Determine noise terms ($\tilde{\sigma}_x$ dBt, $\tilde{\sigma}_z$ dBt) such that the associated diffusion coefficient $\tilde{a}_{zz}$ and Itô-Stokes drift $\partial_z \tilde{a}_{xz}$ match or approximate the given parameterizations ($a_{zz}$, $u_s$).

* Typically, this involves solving optimization problems for weighted coefficients with a predefined set of basis functions

* Here, we consider a particular solution projected onto a set of localized basis functions {$e_n$} in $L^2([-H, 0], \mathbb{R})$ (e.g., regularized Dirac functions, wavelets):

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

$$ \tilde{\sigma}_x(z) dB_t = \sqrt{2} \sum_n \langle a_{zz}^{-1/2} u_s, e_n \rangle e_n(z) d\beta_t^n $$

$u_s = \int_{-H}^2 u_s(\zeta) d\zeta = u_s/(2k)$, $\beta_t^n$: standard Brownian motions (independent to $W_t$)

* Stronger winds correspond to increased vertical noise, while smaller waves enhance horizontal noise at the surface

<figure><img src="image_8.png" /></figure>
1. Introduction

2. Stochastic Formulation

3. Uncertainty Representation

4. Statistical Analyses

5. Conclusion
* Parametrized PDE (p-PDE) [McWilliams et al., 1997]:

$$ \partial_t \mathbf{u}_e = -if(\mathbf{u}_e + \mathbf{u}_s) + \partial_z(a_{zz} \partial_z \mathbf{u}_e) $$

* Parametrized stochastic PDE (p-SPDE):

$$ d\mathbf{u}_e = (\partial_z(a_{zz}\partial_z(\mathbf{u}_e + \mathbf{u}_s)) - if(\mathbf{u}_e + \mathbf{u}_s))dt - (\sigma_z\partial_z(\mathbf{u}_e + \mathbf{u}_s) + if\sigma_x)dB_t $$

* With the same type of boundary condition and initial condition:

$$ \rho_0 a_{zz} \partial_z \mathbf{u}_e|_{z=0} = \tau_w, \quad \rho_0 a_{zz} \partial_z \mathbf{u}_e|_{z=-H} = 0, \quad \mathbf{u}_e(t=0) = 0 $$

<table><thead><tr><th>Parameters</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td><i>f</i></td><td>0.73 × 10<sup>-4</sup> s<sup>-1</sup></td><td>Coriolis freq.</td></tr><tr><td><b>u&#x0302;<sub>a</sub></b></td><td>(5, 0) m s<sup>-1</sup></td><td>Mean wind</td></tr><tr><td>Σ<sub>a</sub></td><td>5 m s<sup>-1</sup></td><td>SD of transient wind</td></tr><tr><td><i>T</i><sub>a</sub></td><td>1 day</td><td>Memory time of wind</td></tr><tr><td>α</td><td>0.8 m</td><td>Surface wave amp.</td></tr><tr><td>λ = 2π/k</td><td>60 m</td><td>Stokes wavelength</td></tr><tr><td>Θ<sub>s</sub></td><td>0°</td><td>Stokes mean angle</td></tr><tr><td>Σ<sub>s</sub></td><td>5°</td><td>SD of Stokes angles</td></tr></tbody></table>
Comparison of Ensemble Statistics

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

* The p-SPDE yields smoother mean profiles, higher uncertainty, increased energy, and stronger correlations than the p-PDE.
Comparison of Higher-Order Statistics

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

* The p-SPDE exhibits more skewed distributions and frequent extreme values than the p-PDE, especially near the surface.
Comparison of evolution after wind stops

**Figure:** The wind is stopped ($\tau_W = 0$) after 30 days of simulation. The mean Ekman spirals at different days in (a, b) are represented by various colors. In each spiral, a darker color indicates a deeper position with a depth step of 5 meters.

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

* The p-SPDE maintains its mean spiral and exhibits higher residual uncertainties than the p-PDE after wind stoppage.
Comparison across transient wind levels

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

* Higher transient winds lead to increased global energy, deeper boundary layers, and stronger Ekman transport with greater angular uncertainty.

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

* Higher wind fluctuations, and thus stronger noise, lead to an increase in extreme events over time, particularly near the surface.
Comparison across Stokes parameters

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

* Increasing wavelength reduces global energy, transport magnitude, and Ekman layer depth, with smaller wavelengths enhancing wave mixing effects.

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

* Global energy and transport magnitude peak when waves align with wind, decreasing symmetrically with rotation; transport angle closely follows wave direction changes.
1. Introduction

2. Stochastic Formulation

3. Uncertainty Representation

4. Statistical Analyses

5. Conclusion
* This novel stochastic formulation is derived to capture the interplay between random motions induced by wind, waves, and currents with enhanced uncertainty representation.

* Numerical results indicate higher kinetic energy, stronger correlations, and more skewed distributions with increased extreme events near the surface compared to traditional models.

* Sensitivity analysis reveals that transient winds and smaller waves increase energy, deepen circulation, and amplify uncertainty, with the strongest effects when waves and wind align.
* Validate the model using reference data from observations or high-quality numerical simulations [McWilliams et al., 2012]

* Analyze wave feedback mechanisms in a 1D fully coupled atmosphere-ocean boundary layer model [Lewis and Belcher, 2004]

* Quantify the influence of wave mixing on the evolution of thermal fronts in a 2D (x, z) stratified Ekman model [Crowe and Taylor, 2018]
# Questions?
<figure><img src="image_16.png" /></figure>

Bauer, W., Chandramouli, P., Chapron, B., Li, L., and Mémin, E. (2020).

Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a
stochastic approach.

J. Phys. Oceanogr., 50(4):983-1003.

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

Crowe, M. N. and Taylor, J. R. (2018).

The evolution of a front in turbulent thermal wind balance. Part 1. theory.

J. Fluid Mech., 850:179-211.

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

Ekman, V. W. (1905).

On the influence of the earth's rotation on ocean-currents.

Ark. Mat. Astr. Fys., 2(11):1-52.

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

Lewis, D. M. and Belcher, S. E. (2004).

Time-dependent, coupled, Ekman boundary layer solutions incorporating Stokes drift.

Dyn. Atmos. Oceans, 37(4):313-351.
<figure><img src="image_20.png" /></figure>

McWilliams, J. C. and Huckle, E. (2006).

Ekman layer rectification.

J. Phys. Oceanogr., 36(8):1646-1659.

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

McWilliams, J. C., Huckle, E., Liang, J.-H., and Sullivan, P. P. (2012).

The wavy Ekman layer: Langmuir circulations, breaking waves, and Reynolds stress.

J. Phys. Oceanogr., 42(11):1793-1816.

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

McWilliams, J. C., Sullivan, P. P., and Moeng, C.-H. (1997).

Langmuir turbulence in the ocean.

J. Fluid Mech., 334:1-30.

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

Phillips, O. M. (1977).

*The Dynamics of the Upper Ocean.*

Cambridge University Press.
<figure><img src="image_24.png" /></figure>

Wenegrat, J. O. and McPhaden, M. J. (2016).

Wind, waves, and fronts: Frictional effects in a generalized Ekman model.

*J. Phys. Oceanogr.*, 46(2):371–394.
* Pseudo-spectral Chebyshev method in space and linear-implicit Euler scheme in time for the SPDE:

$$ \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) $$

$\mathbf{u}_j^n := \mathbf{u}(z_j, t_n)$, $A = -ifI + DaD$, $D$: Chebyshev differentiation matrix, $F$, $G$: RHS operators.

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

* Exact scheme for the Markovian wind process:

$$ \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 $$

$\bar{\mathbf{u}}_a$: mean wind, $T_a$: memory time, $\xi$: standard normal variables.

* Code: https://github.com/matlong/SGESM