# Correction of GLS turbulence closure for wave-resolving models with stratification

Patrick Marchesiello, Simon Treillou

► To cite this version:

Patrick Marchesiello, Simon Treillou. Correction of GLS turbulence closure for wave-resolving models with stratification. Ocean Modelling, 2023, 184, 102212 [7 p.]. 10.1016/j.ocemod.2023.102212. hal-04720383

HAL Id: hal-04720383

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

Submitted on 20 Mar 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
# Correction of GLS turbulence closure for wave-resolving models with stratification

Patrick Marchesiello<sup>a,*</sup>, Simon Treillou<sup>a</sup>

<sup>a</sup>IRD/LEGOS, University of Toulouse, France

## Abstract

We build on the work of Mayer and Madsen (2000) and Larsen and Fuhrman (2018) to correct for overmixing in wave potential flow regions due to the linear instability of two-equation closure systems. We extend it to the Generic Length Scale (GLS) turbulent closure approach, while testing it for the first time on a stratified nearshore problem. The modified turbulence equations affect only the nonbreaking wave potential flow region where the vorticity is less than 10% of the strain rate and thus can preserve stratification on the innershelf under favorable heat and wind forcing conditions. This work is a necessary step before engaging in full 3D wave-resolving studies of surf-shelf exchange.

**Keywords:** Nearshore, Wave-resolving model, Wave-induced mixing, Turbulence closure

## 1. Introduction

Several 3D wave-resolving, free-surface and terrain-following models based on the Reynolds-Averaged Navier-Stokes (RANS) equations have emerged for the nearshore zone, e.g., SWASH (Zijlema et al., 2011), NHWAVE (Derakhti et al., 2016) and CROCO (Marchesiello et al., 2021). In these models, the explicit overturning of the free surface is excluded and the breaking wave is modeled with a single-valued free surface which behaves as a dissipating bore. Despite the absence of explicit overturning (replaced by parameterized turbulence), these models can be accurate as well as computationally efficient in the study of waves and wave-driven mean and transient circulation. They can thus be applied to the entire

<footnote>*Corresponding author
Email address: patrick.marchesiello@ird.fr (Patrick Marchesiello)</footnote>
3D nearshore zone, including the innershelf area and be coupled with biogeochemical or sediment transport models to explicitly simulate the effects of realistic short-crested and asymmetric waves (Marchesiello et al., 2021, 2022). But stratification or Earth rotation has not been included in these earlier works, even though they may be important on the shelf and in the exchange mechanisms with the surfzone.

In previous studies (Marchesiello et al., 2021, 2022), the Coastal and Regional Ocean COmmunity model (CROCO) was used with the *k*-ω or *k*-ε two-equation turbulence closure models as part of a Generic Length Scale (GLS) method (Warner et al., 2005). The *k*-ω model generally performed best in the surfzone without modification of the standard model parameters from Wilcox (1988). This is a conclusion shared by several modeling studies since Mayer and Madsen (2000), possibly due to a more accurate near-wall treatment by *k*-ω models (Devolder et al., 2018), which is needed in the wall-bounded turbulent shear flows of the surfzone (Trowbridge and Elgar, 2001). However, both models suffer from overmixing in areas of potential (irrotational) flow, typically in regions outside the surfzone where waves are not breaking. This problem has long been identified (Lin and Liu, 1998; Mayer and Madsen, 2000; Bradford, 2011; Brown et al., 2016; Devolder et al., 2018; Larsen and Fuhrman, 2018) and understood as a kind of paradox where vorticity, rather than strain rate, should be the source of turbulence, even though the role of strain rate is theoretically justified by the Reynolds averaging process (Kato and Launder, 1993). The generation of turbulence in unexpected regions may thus appear to challenge the validity of Boussinesq eddy viscosity hypothesis on the local relationship between mean and turbulent quantities (Schmitt, 2007), which led in particular to the realizability constraints to the introduction of a realizability constraint based on physical assumptions (Durbin, 2009).

Recently, however, Larsen and Fuhrman (2018), building on earlier analysis of Mayer and Madsen (2000) that went rather unnoticed, have better identified the causes of excessive vertical mixing in areas of potential flow, revealing a linear instability of the two-equation closure system. It is the interaction between the two turbulent prognostic variables in the system that can exponentially amplify the generation of turbulent energy in unexpected regions. The resulting overmixing can affect the amplitude of waves propagating to the
39 nearshore zone and degrade the prediction of the wave breaking point (Lin and Liu, 1998; 40 Bradford, 2000; Mayer and Madsen, 2000). We found in the present study that, perhaps more 41 importantly, overmixing can also affect the maintenance of stratification on the innershelf. 42 The last point is generally overlooked because 3D wave-resolving models have been primarily 43 focused on the surfzone without considering exchanges with the innershelf, in contrast to 44 recently emerging wave-averaged model studies (Kumar and Feddersen, 2017; Wang et al., 45 2021), which do not suffer from the overmixing problem. With advances in non-hydrostatic 46 models and computational resources, the study of surf-shelf exchange — describing the fate 47 of pollutants, biological matter, and sediments in the nearshore zone (e.g., Brown et al. 48 2015) — can now also be accessed by representing phase-resolved waves and stratification. 49 However, to proceed in this direction, the problem of excessive mixing in nonbreaking waves 50 must first be corrected.

51 Note that in wave-averaged models, a certain amount of nonbreaking wave-induced mix-
ing is actually introduced by semi-empirical formulas (Qiao et al., 2004; Ghantous and
52 Babanin, 2014). It is much smaller than that caused by wave breaking but it reaches deeper
53 and can affect the sea surface temperature of ocean general circulation models (Babanin,
54 2023). The wave-induced turbulence involved in these parameterizations assumes various
55 processes, such as interaction with pre-existing turbulence. Here, the generation of intense
56 innershelf turbulent mixing appears to be more of an artifact of the turbulent models under
57 conditions of nearly potential flow. To address this specific issue, we build on the works of
58 Mayer and Madsen (2000) and Larsen and Fuhrman (2018) and extend it to a more generic
59 approach within GLS, while testing for the first time the effect that a stabilized turbulence
60 model has on innershelf stratification in an idealized nearshore configuration.
61

## 2. Modification of the turbulence closure model

62 CROCO is a nonhydrostatic, free-surface and terrain-following model, developed around 63 the Regional Oceanic Modeling System (Shchepetkin and McWilliams, 2005; Debreu et al., 64 2012). Its capabilities include high-performance computation of high-order-accurate dis- 65 cretized equations, and coupling with atmosphere, waves, biogeochemistry, sediment and 66
turbulence models. It has been applied to a variety of configurations, from regional and shelf circulations to very fine-scale processes, such as wave-induced nearshore circulation (Marchesiello et al., 2015, 2021, 2022). It thus appears naturally suited for addressing surf-shelf exchange processes in a 3D, rotating and stratified framework. However, innershelf stratification cannot be preserved if the turbulence closure equations produce excessive mixing, as is often the case for RANS wave-resolving models (Larsen and Fuhrman, 2018).

## 2.1. Turbulence closure

The treatment of breaking waves in CROCO involves both a shock-capturing advection scheme (WENO5) and a second-order turbulence closure system — a $k-\epsilon$ or $k-\omega$ model solving the equations for turbulent kinetic energy $k$ and dissipation $\epsilon$ or dissipation rate $\omega \propto \epsilon k^{-1}$ — as part of a Generic Length Scale (GLS) method (Warner et al., 2005). In the absence of buoyancy forcing, the turbulence equations express a balance between transport, diffusion, shear production and dissipation:

$$ \frac{\partial \rho k}{\partial t} = -\vec{\nabla} \cdot (\rho \vec{v} k) + \mathcal{D}_k + \rho (P - \epsilon) \quad (1) $$

$$ \frac{\partial \rho \epsilon}{\partial t} = -\vec{\nabla} \cdot (\rho \vec{v} \epsilon) + \mathcal{D}_{\epsilon} + \rho \frac{\epsilon}{k} (c_{\epsilon 1} P - c_{\epsilon 2} \epsilon) \quad (2) $$

$$ \text{or} \quad \frac{\partial \rho \omega}{\partial t} = -\vec{\nabla} \cdot (\rho \vec{v} \omega) + \mathcal{D}_{\omega} + \rho \frac{\omega}{k} (c_{\omega 1} P - c_{\omega 2} \epsilon) \quad (3) $$

The eddy viscosity $\nu_t = c_\mu l k^{\frac{1}{2}}$ is derived from these equations, with coefficient $c_\mu$ dependent on stability functions, and mixing length $l \propto k^{\frac{3}{2}}\epsilon^{-1}$. According to the Boussinesq eddy viscosity hypothesis, the shear production term for $k$ is $P = 2\nu_t S^2$, with $S^2 = S_{ij}S_{ij}$ and $S_{ij} = \frac{1}{2}(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})$ (using Einstein notation), which is the mean 3D strain rate tensor<fnref n="1" />. All turbulence model parameters are given in Warner et al. (2005), based on Burchard et al. (1998) for $k-\epsilon$ and Wilcox (1988) for $k-\omega$. The turbulence models perform very well in the surfzone, compared with flume experiments (Marchesiello et al., 2021, 2022). However, as in previous studies, they tend to produce excessive mixing in potential flow regions, i.e., on the innershelf.

<footnote n="1">¹The use of a 3D strain rate tensor differs from the GLS approach presented in Warner et al. (2005), where energy production is only due to vertical shear of horizontal flows.</footnote>
## 2.2. Correction in potential flow regions

In many studies of breaking waves (see the literature reviews by Devolder et al. 2018; Larsen and Fuhrman 2018), the predicted turbulence levels from two-equation closures are much higher than what has actually been measured. While many attempts were made for limiting mixing in supposed low-turbulence regions, Larsen and Fuhrman (2018) more specifically identify the causes of overmixing in areas of potential (irrotational) flow, typical of nonbreaking waves. The problem is common to the entire class of 2-equation closure models ($k-\omega$, $k-\epsilon$, and variants) because it is related to an unconditional linear instability of the equation system under low vorticity conditions.

For the $k-\omega$ model, a linearized homogeneous system suitable for the potential flow region beneath surface waves (where a fixed strain rate $S$ and negligible advection and diffusion can be assumed) is:

$$ \frac{\partial k}{\partial t} = P - \epsilon \qquad (4) $$

$$ \frac{\partial \omega}{\partial t} = \frac{\omega}{k}(c_{\omega 1}P - c_{\omega 2}\epsilon) \qquad (5) $$

Noting that the specific dissipation rate $\omega$ is related to the energy dissipation $\epsilon$ by $\epsilon = c_\mu k \omega$ and that the eddy viscosity $\nu_t = \frac{k}{\omega}$, we can rewrite the system of equations in terms of only the variables $k$ and $\omega$:

$$ \frac{\partial k}{\partial t} = 2\frac{k}{\omega}S^2 - c_{\mu}k\omega \qquad (6) $$

$$ \frac{\partial \omega}{\partial t} = 2c_{\omega 1}S^2 - c_{\omega 2}c_{\mu}\omega^2 \qquad (7) $$

Considering as Larsen and Fuhrman (2018) that $\omega$ tends asymptotically to $\omega_\infty$ as $\partial\omega/\partial t$ tends to 0 and replacing $\omega$ by $\omega_\infty$ in the $k$ equation gives an exponential growth $k \sim e^{\sigma_\infty t}$ with $\sigma_\infty$ a fraction of $S$ always positive, i.e., the model is unconditionally unstable.

The instability of this system thus corresponds to an exponential increase in the turbulent energy $k$, particularly in areas of potential flow where dissipation is never strong enough to control it. Larsen and Fuhrman (2018) propose to increase dissipation curb the excessive energy production where the strain rate $S$ is much greater than the vorticity $\Omega$ (defined
by $\Omega^2 = \Omega_{ij}\Omega_{ij}$ with $\Omega_{ij} = \frac{1}{2}(\frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i})$, the rotation rate tensor) and demonstrate that this correction stabilizes the solution of the linear system as $\sigma_\infty$ becomes negative. The limitation is written as a majoration of $\omega$ in the viscosity $\nu_t = \frac{k}{\omega^*}$ used in the production term $P = 2\nu_t S^2$, with the following function:

$$ \omega^* = \max \left( 1, \lambda \frac{c_{\omega 2} S^2}{c_{\omega 1} \Omega^2} \right) \omega \qquad (8) $$

where $\lambda$ is a limitation parameter (taken as 0.05 in Larsen and Fuhrman 2018). A similar condition is given for the $k-\epsilon$ model ($\nu_t = c_\mu \frac{k^2}{\epsilon^*}$), where:

$$ \epsilon^* = \max \left( 1, \lambda \frac{c_{\epsilon 2} S^2}{c_{\epsilon 1} \Omega^2} \right) \epsilon \qquad (9) $$

This correction is reminiscent of the one proposed for the stagnation point anomaly by Kato and Launder (1993), in a different study context, where the turbulent energy production $P = 2\nu_t S^2$ is replaced by $P = 2\nu_t S\Omega$. It also resembles that of Mayer and Madsen (2000) who take $P = 2\nu_t \Omega^2$, thus making production vanish in pure straining flow. According to Larsen and Fuhrman (2018), these attempts significantly reduce the growth of instability in the potential flow region, but do not resolve its formal unconditionality. In addition, in our numerical experiments, the use of vorticity to replace all or part of the strain rate in the production term $P$ tends to under-represent the production of turbulent energy in the surfzone, which is undesirable. We thus concur with Mayer and Madsen (2000) and Durbin (2009) that altering the strain rate tensor in $P$ is not a satisfactory solution and we choose instead to majorate $\epsilon$ (when computing $\nu_t$) to control the instability of the turbulence equation system. The latter solution may appear merely as an alternative to limiting the mixing length $l$ as in Bradford (2011) or the coefficient $c_\mu$ as in realizable models (Shih et al., 1995), but a limitation based on the ratio $S/\Omega$ makes it a more effective constraint on the instability (Fuhrman and Li, 2020).

## 2.3. Implementation in GLS

For the implementation in the CROCO GLS code, we impose the Larsen and Fuhrman (2018) limitation on the variable $\epsilon$ whatever the choice of closure, since the turbulent dissipation $\epsilon$ is computed in all types of two-equation models handled by GLS. Starting from
the k-ω model, we use the relation between $\epsilon$ and $\omega$ ($\epsilon = c_\mu k\omega$) to assume a limitation of the form:

$$ \epsilon^* = \max \left( 1, \lambda \frac{c_{\omega 2} S^2}{c_{\omega 1} \Omega^2} \right) \epsilon \quad (10) $$

This expression is similar to Eq. 9 for the $k-\epsilon$, apart from the factor $c_{\omega 2}/c_{\omega 1}$ replacing $c_{\epsilon 2}/c_{\epsilon 1}$. However, the fraction $c_2/c_1$ does not vary much with the choice of closure model (Warner et al., 2005) and is 1.2-1.5. Considering the uncertainty in choosing $\lambda$, we can propose a generic limitation:

$$ \epsilon^* = \max \left( 1, \lambda \frac{S^2}{\Omega^2} \right) \epsilon \quad (11) $$

In our numerical experiments (next section), this formulation controlled the growth of instability on the innershelf with $\lambda \sim 0.1$. In this case, the correction is triggered as soon as the vorticity is less than about 10% of the strain rate. We also note that increasing the power law of the correction function stiffens the intensity of the correction below the 10% threshold and further decreases the mixing on the innershelf. Therefore, in the following, we will test the modified formulation:

$$ \epsilon^* = \max \left( 1, \lambda \frac{S^4}{\Omega^4} \right) \epsilon \quad (12) $$

with a baseline $\lambda$ value of 0.01.

# 3. Numerical experiments

## 3.1. Ideal test case

The idealized model configuration is representative of a nearly planar beach representative of southern California (Kumar and Feddersen, 2017) or elsewhere with a smooth sandbar that triggers wave breaking, and a cross-shore circulation (Fig. 1). Details of the numerical implementation of CROCO for wave-resolving simulations is given in (Marchesiello et al., 2021, 2022), only the configuration is presented here. For simplicity, we use a 2D (x-z) channel with no alongshore dimension. Depths range from 7 m offshore to 1 m above mean sea level at the coast and a gentle sandbar is present about 100 m from shore. The shoreline position is allowed to vary with swash oscillation, based on a wetting-drying scheme (Warner
et al., 2013). Initially, the water is at rest and the stratification is given by a temperature profile $T(z) = 8 + 10e^{z/50}$ (Fig. 2), roughly representing the mid-latitude summer situation.

A monochromatic wave is generated at the offshore boundary with 0.5 m wave amplitude and 10 s period. A non-slip condition is imposed on the boundaries of the channel sidewalls. The horizontal grid spacing is 0.5 m, the vertical grid has 10 vertical levels refined at the surface and bottom and the model time step is 0.025 s. For bed shear stress, the logarithmic law of the wall is used with a roughness length of 1 mm.

The comparison of velocity $u$, temperature $T$ and eddy viscosity $\nu_t$ at 30 min into after the start of the simulation is presented in Figure 2. It can be seen that without limitation, the eddy viscosity produced by the $k-\omega$ model increases to values well above 0.01 m$^2$/s and the high viscosity signal propagates from the surfzone to the shelf and gradually invade the entire model area (Fig. 3). With the corrected closure system, the viscosity does not develop outside the surfzone, i.e., in the potential flow region where waves do not break. In this case, the viscosity in the surfzone has the expected structure and magnitude ($\nu_t = 0.01-0.05$ m$^2$/s; Battjes 1975). It is important to note that in the breaker zone, the viscosity — and related velocity profile — is similar in both the standard and corrected cases, showing that the limitation does not seem to affect the assumed turbulent zone much (which is desirable), but this will be verified against laboratory data in the next section.

The effect of vertical mixing on stratification is dramatic (Fig. 2). At viscosity values of about 0.01 m$^2$/s and above, mixing can break up the stratification in a matter of minutes ($\tau_{mix} = l^2/\nu_t$ with $l$ limited by bottom depth). With the original $k-\omega$ model, such high viscosity values reach the entire domain within half an hour and stratification breakdown follows everywhere. On the contrary, stratification is maintained on the innershelf outside the surfzone with the corrected turbulence model. In this case, the result becomes similar to those presented by wave-averaged models (lacking a process to generate a nonbreaking wave eddy viscosity; Kumar and Feddersen 2017) and to nearshore observations (Sinnett and Feddersen, 2019), where spring-summer stratification is preserved outside the surfzone as long as winds are sufficiently light.

Next, we checked that the corrections made on the turbulent dissipation $\epsilon$ are working
not only for the *k*-ω model but also for the *k*-ε and the generic model of Umlauf et al. (2003) which is referred to as *gen* model in Warner et al. (2005). Figure 4 shows the result of the three models. It confirms the report of Marchesiello et al. (2021) that the *k*-ε model produces less turbulent energy than the *k*-ω model, with as a consequence an over-prediction of velocity shear and undertow magnitude. The *gen* model gives intermediate results. Most importantly here, the correction has worked in all cases and no exponential growth of turbulent energy and viscosity has occurred outside of the surfzone. Our correction of turbulent dissipation thus appears appropriate for a GLS approach.

## 3.2. Surfzone verification with LIP experiment

Validation of surfzone currents by comparison with flume experiments are provided in Marchesiello et al. (2021, 2022). To make sure that the change of mixing due to the GLS corrections does not affect the surfzone excessively, we present here again some comparison with the full-scale LIP experiment of Roelvink and Reniers (1995). Following Durbin (2009), we assume that since eddy viscosity models are intended for prediction of mean flow, mean flow data are a good basis for calibration.

In the LIP1B experiment, a JONSWAP wave spectrum with significant wave height $H_s = 1.4$ m and period $T_p = 5$ is generated that produces wave breaking at various position around a sandbar. The horizontal grid resolution is 25 cm with 20 vertical levels refined near the bottom. For more details on the LIP experiment and model configuration, the reader is referred to the previous reports cited above.

Figure 5 shows a comparison of the model cross-shore currents with data, and the match appears very good throughout the complex morphology of the beach. The waves break everywhere depending on their height but the breaking is more intense on the sandbar, where the undertow has a strong shear and maximum intensity of about 30 cm/s. The mean absolute error of all experiments (Fig. 5, bottom panel) is 3-4 cm/s (RMSE=2-3 cm), close to the measurement error of 2 cm/s. As discussed in Marchesiello et al. (2021), these results show the validity of a RANS approach for estimating breaking-induced turbulence (via parameterization) and calculating lower frequency circulation explicitly.
The experiments with correction of the $k-\omega$ model do not modify much the results, but differences appear between the case with strong limitation ($\lambda = 0.1$) and moderate limitation ($\lambda = 0.01$). In the former, the mean error with measurements is 4.9 cm while the latter is 3.7 cm, close to 3.6 cm for the standard case. The error in the strong limitation case is most significant in the offshore region of the sandbar where only the largest waves in the spectrum break, implying intense but transient turbulent energy. However, the error reduces again at the most offshore location where turbulent energy is supposedly weakening. The moderate limitation case performs well and even reduces the error, relative to the standard model, at a few inshore locations as well as at the most offshore location. The $\lambda$ value of 0.01 confirms to be a good choice for an optimal correction level. By further decreasing $\lambda$, the instability reappears.

The genericity of the $\lambda$ value is unclear. We have increased the power law of the limiting function (Sec. 2.3) in order to make the correction more robust and less sensitive to the limiting value. The fact that the correction works for different GLS model types seems to confirm this idea. We believe that the method is robust, at least for the problem of waves propagating in the nearshore zone. We also verified that adding wind stress to the model (in the idealized nearshore configuration of Sec. 3.1) allows mixing of the stratified innershelf as expected, i.e. the correction does not prevent the generation of mixing outside the surfzone if the forcing conditions are favorable.

# 4. Conclusion

In this paper, we build on the work of Mayer and Madsen (2000) and Larsen and Fuhrman (2018) to correct for overmixing in potential flow regions due to the linear instability of two-equation closure systems. We extend it to a more generic GLS approach, while testing for the first time the corrected closure model on a stratified nearshore problem. The corrected GLS formulation changes little in the surfzone where high vorticity and turbulence are expected, while successfully addressing the problem of overmixing in regions where the waves are not breaking, allowing the model to maintain stratification on the shelf provided the vorticity is less than 10% of the strain rate. This work is a necessary step before engaging in full wave-
232 resolving 3D studies of surf-shelf exchange that involve transient rip currents, stratification
233 and also Earth rotation. This will be possible due to advances in non-hydrostatic free-surface
234 models, which are emerging as a cost-effective and accurate alternative to models that fully
235 resolve the wave breaking process (Marchesiello et al., 2021).

236 Acknowledgement

237 This research has received support from a consortium of French research agencies, as part
238 of CROCO's development project (Insu GdR n°2014 named CROCO) and from the French
239 Naval Hydrographic and Oceanographic Service (SHOM DGA-Protevs II). We thank Dano
240 Roelvink for sharing the LIP data. Apart from these, all data were acquired by the authors
241 and the CROCO source code is freely available at www.croco-ocean.org. Both laboratory
242 and modeling data are available upon request.

243 References

244 Babanin, A.V., 2023. Ocean waves in large-scale air-sea weather and climate systems. Journal of Geophysical
245 Research: Oceans 128, e2023JC019633.

246 Battjes, J., 1975. Modeling of turbulence in the surfzone, in: A.S.C.E., New York, N. (Ed.), Symposium on
247 Modeling Techniques, pp. 1050–1061.

248 Bradford, S., 2000. Numerical simulation of surf zone dynamics. Journal of Waterway, Port, Coastal, and
249 Ocean Engineering 126, 1–13.

250 Bradford, S., 2011. Nonhydrostatic model for surf zone simulation. Journal of Waterway, Port, Coastal, and
251 Ocean Engineering 137, 163–174.

252 Brown, J.A., MacMahan, J.H., Reniers, A.J.H.M., Thornton, E.B., 2015. Field observations of surf
253 zone–inner shelf exchange on a rip-channeled beach. Journal of Physical Oceanography 45, 2339 – 2355.

254 Brown, S., Greaves, D.M., Magar, V., Conley, D.C., 2016. Evaluation of turbulence closure models under
255 spilling and plunging breakers in the surf zone. Coastal Engineering 114, 177–193.

256 Burchard, H., Petersen, O., Rippeth, T.P., 1998. Comparing the performance of the mellor-yamada and the
257 $k-\epsilon$ two-equation turbulence models. Journal of Geophysical Research: Oceans 103, 10543–10554.

258 Debreu, L., Marchesiello, P., Penven, P., Cambon, G., 2012. Two-way nesting in split-explicit ocean models:  
Algorithms, implementation and validation. Ocean Modelling 49–50, 1 – 21.

260 Derakhti, M., Kirby, J.T., Shi, F., Ma, G., 2016. Nhwave: Consistent boundary conditions and turbulence modeling. Ocean Modelling 106, 121 – 130.
261
Devolder, B., Troch, P., Rauwoens, P., 2018. Performance of a buoyancy-modified $k-\omega$ and $k-\omega$ sst turbulence model for simulating wave breaking under regular waves using openfoam. Coastal Engineering 138, 49–65.

Durbin, P.A., 2009. Limiters and wall treatments in applied turbulence modeling. Fluid Dynamics Research 41, 012203.

Fuhrman, D.R., Li, Y., 2020. Instability of the realizable $k-\epsilon$ turbulence model beneath surface waves. Physics of Fluids 32, 115108.

Ghantous, M., Babanin, A.V., 2014. One-dimensional modelling of upper ocean mixing by turbulence due to wave orbital motion. Nonlinear Processes in Geophysics 21, 325–338.

Kato, M., Launder, B., 1993. The modelling of turbulent flow around stationary and vibrating square cylinders, pp. 1–6.

Kumar, N., Feddersen, F., 2017. The effect of stokes drift and transient rip currents on the inner shelf. part i: No stratification. Journal of Physical Oceanography 47, 227–241.

Larsen, B.E., Fuhrman, D.R., 2018. On the over-production of turbulence beneath surface waves in reynolds-averaged navier–stokes models. Journal of Fluid Mechanics 853, 419–460.

Lin, P., Liu, P.L.F., 1998. A numerical study of breaking waves in the surf zone. Journal of Fluid Mechanics 359, 239–264.

Marchesiello, P., Auclair, F., Debreu, L., McWilliams, J., Almar, R., Benshila, R., Dumas, F., 2021. Tridimensional nonhydrostatic transient rip currents in a wave-resolving model. Ocean Modelling 163, 101816.

Marchesiello, P., Benshila, R., Almar, R., Uchiyama, Y., McWilliams, J.C., Shchepetkin, A., 2015. On tridimensional rip current modeling. Ocean Modelling 96, 36–48. Waves and coastal, regional and global processes.

Marchesiello, P., Chauchat, J., Shafiei, H., Almar, R., Benshila, R., Dumas, F., Debreu, L., 2022. 3d wave-resolving simulation of sandbar migration. Ocean Modelling 180, 102127.

Mayer, S., Madsen, P.A., 2000. Simulation of Breaking Waves in the Surf Zone using a Navier-Stokes Solver. pp. 928–941.

Qiao, F., Yuan, Y., Yang, Y., Zheng, Q., Xia, C., Ma, J., 2004. Wave-induced mixing in the upper ocean: Distribution and application to a global ocean circulation model. Geophysical Research Letters 31.

Roelvink, J.A., Reniers, 1995. IP 11D delta flume experiments : a dataset for profile model validation. WL / Delft Hydraulics.

Schmitt, F.G., 2007. About boussinesq’s turbulent viscosity hypothesis: historical remarks and a direct evaluation of its validity. Comptes Rendus Mécanique 335, 617–627. Joseph Boussinesq, a Scientist of bygone days and present times.

Shchepetkin, A.F., McWilliams, J.C., 2005. The regional oceanic modeling system (roms): a split-explicit,
296 free-surface, topography-following-coordinate oceanic model. Ocean Modelling 9, 347-404.

297 Shih, T.H., Liou, W.W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new k-ε eddy viscosity model for high reynolds number turbulent flows. Computers and Fluids 24, 227-238.

299 Sinnett, G., Feddersen, F., 2019. The nearshore heat budget: Effects of stratification and surfzone dynamics. Journal of Geophysical Research: Oceans 124, 8219-8240.

301 Trowbridge, J., Elgar, S., 2001. Turbulence measurements in the surf zone. Journal of Physical Oceanography 31.

303 Umlauf, L., Burchard, H., Hutter, K., 2003. Extending the k-ω turbulence model towards oceanic applications. Ocean Modelling 5, 195-218.

305 Wang, P., McWilliams, J.C., Uchiyama, Y., 2021. A nearshore oceanic front induced by wave streaming. Journal of Physical Oceanography 51, 1967-1984.

307 Warner, J., Sherwood, C., Arango, H., Signell, R., 2005. Performance of four turbulence closure methods implemented using a generic length scale method. Ocean Modelling 8, 81-113.

309 Warner, J.C., Defne, Z., Haas, K., Arango, H.G., 2013. A wetting and drying scheme for ROMS. Computers and Geosciences 58, 54-61.

311 Wilcox, D.C., 1988. Reassessment of the scale-determining equation for advanced turbulence models. AIAA Journal 26, 1299-1310.

313 Zijlema, M., Stelling, G., Smit, P., 2011. Swash: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters. Coastal Engineering 58, 992-1012.
<figure><img src="image_2.png" /><figcaption>Figure 1: Snapshot of CROCO wave height and cross-shore velocity $u$ propagating across an idealized nearshore configuration. The simulation is done at 50 cm resolution and with the $k-\omega$ turbulence closure model.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Wave-mean velocity $u$, temperature $T$ and eddy viscosity $\nu_t$ at initial time (left) and 30 min into the simulation (center and right panels). Center: solution of the standard $k-\omega$ closure model; right: solution of the corrected closure model.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Evolution of eddy viscosity $\nu_t$ 10, 20 and 30 min after the start of the simulation. Left: solution resulting from the standard k-ω closure model; right: solution from the corrected closure model.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Eddy viscosity $\nu_t$ 30 min into the simulation from three different versions of the GLS closure model: k-ω, k-ε and gen.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Model simulation of the large-scale LIP11-1B Flume experiment. Top) Comparison of simulated and measured current profiles $u$ and sensitivity to the k-ω turbulence model correction ($\lambda = 0.1$ and $0.01$); Bottom) Mean Absolute Error of $u$ relative to measurements for all experiments.</figcaption></figure>