Non-breaking Waves Research Articles

Impulse Wave Runup on Steep to Vertical Slopes

Impulse waves are generated by landslides or avalanches impacting oceans, lakes or reservoirs, for example. Non-breaking impulse wave runup on slope angles ranging from 10° to 90° (V/H: 1/5.7 to 1/0) is investigated. The prediction of runup heights induced by these waves is an important parameter for hazard assessment and mitigation. An experimental dataset containing 359 runup heights by impulse and solitary waves is compiled from several published sources. Existing equations, both empirical and analytical, are then applied to this dataset to assess their prediction quality on an extended parameter range. Based on this analysis, a new prediction equation is proposed. The main findings are: (1) solitary waves are a suitable proxy for modelling impulse wave runup; (2) commonly applied equations from the literature may underestimate the runup height of small wave amplitudes; (3) the proposed semi-empirical equations predict the overall dataset within ±20% scatter for relative wave crest amplitudes ε, i.e., the wave crest amplitude normalised with the stillwater depth, between 0.007 and 0.69.

Towards a more complete tool for coastal engineering: solitary wave generation, propagation and breaking in an SPH-based model

ABSTRACTThe present work describes the implementation of an advanced solitary wave generation system in the mesh-less SPH-based DualSPHysics model to simulate tsunami-like solitary waves. Three different generation theories have been implemented and are extended to generate multiple solitary waves. The numerical model is validated against theoretical solutions and physical model results from three different experimental campaigns, which investigated respectively: i) the forces exerted on harbour protections; ii) the run-up of solitary waves on a gentle beach; and iii) the impact of double solitary waves on two cylindrical reservoirs. The differences of modelling breaking and non-breaking wave conditions are highlighted. The results demonstrate the capability of DualSPHysics to represents the main hydrodynamic properties of solitary waves when interacting with the shoreline or with coastal structures. The applicability of DualSPHysics to coastal engineering problems is enhanced, being the model already capable of generating monochromatic waves and random sea states.

Numerical Evaluation of Design Rules for Non-Breaking Wave Loads on Vertical Walls

This paper describes a numerical evaluation of design rules for the determination of wave loads of non-breaking waves on vertical structures. Design guidelines were proposed by Sainflou (1928) and Goda (1974). These guidelines use geometric parameters of the structure, an incident wave height and a wave period. In practice (cf. CERC, 1984), a Rayleigh distribution of individual wave heights is assumed to determine the design wave height in an irregular wave field. Their reliability and range of applicability are poorly known, especially when the incident wave condition consists of a mixed sea state, like a local wind sea and a low-frequency (swell) component. To validate the above described design methods, we applied the non-hydrostatic numerical wave model SWASH to simulate wave loading on a rigid vertical wall for single and mixed sea states. In addition, we compared the results with linear wave theory and the spectral response approach using transfer functions based on linear wave theory.

LONG-TERM COASTAL EVOLUTION MODELLING OF LONGSHORE BARS

An extended version of a numerical model introduced by Larson et al. (2013) to simulate long-term cross shore material exchange for the subaqueous portion of the profile has been developed. Efforts have focused on improving the model to better account for beach systems consisting of two bars (inner and outer bar), as well as simulating the feeder response over time of nearshore dredged material bars, intended to function as beach nourishment. The theory for the evolution of a single-bar to a two-bar system was modeled, considering an inner and an outer bar, where the outer bar is of primary interest with the purpose of predicting the behavior of placed dredged material. The cross-shore sediment transport rate is based on the evolution equation for the bar system response to the hydrodynamic forcing by reference to its equilibrium condition, where the change in the bar volume is based on a set of wave criteria, describing the onset of a new breaking zone when the outer bar forms. Empirical formulas are employed for the bar equilibrium volume and for coefficients determining the bar response rate. In this study, a description of the extended model and the results from the model component validation at two different sites in USA (Duck, North Carolina, and Cocoa Beach, Florida) are presented. Duck measurements have detected that some bars form in the nearshore and move all the way offshore (eventually deflating by non-breaking waves). At the same time, it was equally observed that a lot of inner bars formed in shallow water do not move offshore but remain as inner bars all the time. According to this, the developed model considers that the inner bar will not become the outer bar, but material previously dedicated to the inner bar will be available for the outer bar. Overall, the present study demonstrates the potential for using rather simple models, based on the definition of an equilibrium state that is compared to the current state and the magnitude of offshore wave forcing to drive the changes in the profile. The methodology employed here allowed to quantitatively reproduce the main trends in the subaqueous beach profile response in a long-term perspective as a function of the bar volumes disequilibrium, the magnitude of the incident wave height and the dimensionless fall velocity to move the sand with a time varying forcing.

WAVE LOADING FOR RECURVED PARAPET WALLS IN NON-BREAKING WAVE CONDITIONS: ANALYSIS OF THE INDUCED IMPULSIVE FORCES

This paper describes 2-D numerical simulations aiming to reproduce the pressure impulse named confined-crest impact (Castellino et al., 2018), which occurs when a recurved parapet wall and non-breaking wave conditions are interacting. The simulations are carried out by using the IH2VOF and IHFOAM, the latter developed as OpenFOAM additional library. The results show a large increase of the pressures and forces value when the recurved part of the vertical parapet results completely occluded by the non-breaking wave crest. A sensitivity analysis has been carried out to study the influence of the geometrical parameters (radius r and opening angle a). It has been found a low variability with respect to the radius increase (from 1.0 m to 2.0 m) and a higher influence related to the opening angle variation. Finally, the non-dimensional force component has been represented as a function of the hydraulic and geometrical parameters by means of the dimensionless product (l/h)*s. These parameters represent the overhang extension seaward of the parapet, the water depth and the wave steepness with reference to deep-water conditions.

HIGHLY-RESOLVED NUMERICAL AND LABORATORY ANALYSIS FOR NONBREAKING SOLITARY WAVE SWASH OVER A STEEP SLOPE

In this paper we study the swash processes generated by a nonbreaking solitary wave running up and down a steep slope (1:3). We use experimental data to study flow features and velocities inside the boundary layer, and numerical modelling to investigate variables not measured during the laboratory experiments, such as pressures and bottom shear stress. We focus on the mechanisms that produce flow separation and vortex formation. Particularly, we study a system of vortices generated under a hydraulic jump during the rundown phase, which was first observed by Matsunaga & Honji (1980).

REDUCTION OF THE WAVE OVERTOPPING DISCHARGE AT DIKES IN PRESENCE OF CROWN WALLS WITH BULLNOSES

New numerical and laboratory investigations on wave overtopping at dikes with crown walls were carried out. The main objective of the experiments, presented for the first time in this contribution, is to investigate the effects of the inclusion of bullnoses on the top of crown walls to reduce the average overtopping discharge q. The study extends the experience available on structures with bullnoses, which is so far limited to dikes with promenades under non-breaking wave conditions. The new data on q resulting from the campaign of experiments are compared with the existing predicting formulae for q of the EurOtop manual (2016), in order to verify and upgrade their range of validity. A formulation for a new correction coefficient γ** to be included in the formulae is proposed to account for the effects of the bullnose also in case of structures subjected to breaking waves. A simple solution to represent the geometry of the bullnoses in the EurOtop Artificial Neural Network (ANN) is investigated. The solution, which avoids the ANN re-training and does not require the inclusion of new input parameters, applied to new and existing data gives promising results.

ROCK ARMOR DAMAGE IN DEPTH-LIMITED BREAKING WAVE CONDITIONS

The armor layer of a mound breakwaters is usually designed with a formula derived from physical tests in non-breaking wave conditions; however, most rubble mound breakwaters are placed in the wave breaking zone where the highest waves break before reaching the structure. The hydraulic stability formulas developed for rock-armored breakwaters in non-breaking conditions are not completely valid to characterize the hydraulic stability of these structures under depth-limited wave attack. In this study, five series of 2D physical tests were carried out on a bottom slope m=1/50 to analyze the hydraulic stability of double-layer rock armored breakwaters in depth-limited breaking wave conditions. Measurements taken by 12 wave gauges placed along the wave flume were compared with estimations of Hm0, H2% and H1/10 obtained from numerical model SwanOne. The significant wave height, Hm0, estimated or measured at a distance 3hs from the toe of the structure was the best characteristic wave to relate armor damage with stability number. The six-power relationship between dimensionless armor damage and stability number, found in this study, explained more than 94% of the variance in the damage observations. This relationship is valid for conventional non-overtopping double-layer rock-armored breakwaters on bottom slope m=1/50 and depth-limited breaking wave conditions.

Cubipod® Armor Design in Depth-Limited Regular Wave-Breaking Conditions

Armor stability formulas for mound breakwaters are commonly based on 2D small-scale physical tests conducted in non-overtopping and non-breaking conditions. However, most of the breakwaters built around the world are located in breaking or partially-breaking wave conditions, where they must withstand design storms having some percentage of large waves breaking before they reach the structure. In these cases, the design formulas for non-breaking wave conditions are not fully valid. This paper describes the specific 2D physical model tests carried out to analyze the trunk hydraulic stability of single- and double-layer Cubipod® armors in depth-limited regular wave breaking and non-overtopping conditions with horizontal foreshore (m = 0) and armor slope (α) with cotα = 1.5. An experimental methodology was established to ensure that 100 waves attacked the armor layer with the most damaging combination of wave height (H) and wave period (T) for the given water depth (hs). Finally, for a given water depth, empirical formulas were obtained to estimate the Cubipod® size which made the armor stable regardless of the deep-water wave storm.

Viscous Flow Past a NACA0012 Foil Below a Free Surface Through the Delta-Plus-SPH Method

The present work is dedicated to the modeling of viscous flow past a NACA0012 foil fixed in a current below a free surface. To this end, the [Formula: see text]-smoothed-particle hydrodynamics (SPH) model has been adopted. This SPH model prevents the inception of the numerical tensile instability in the flow region characterized by negative pressure since a tensile instability control (TIC) has been included. In the TIC, a pressure differencing formulation (PDF) has been adopted for the momentum equation in the flow region characterized by negative pressure. In order to completely remove the numerical noise in the vorticity field, in this work, the PDF is also applied for the region with positive pressure, but except for the free-surface region in order to ensure the free surface stability when wave breaking occurs. The mechanism of PDF being able to eliminate the numerical noise in the vorticity field is also briefly analyzed. In order to reduce the nonconservation of total momentum induced by the PDF, a particle-shifting technique (PST) is implemented in each time step for regularizing the particle position. In the numerical results, [Formula: see text]-SPH results are validated by the experimental data and other verified numerical results. Improvements of the results of [Formula: see text]-SPH with PDF with respect to the ones without using PDF are demonstrated. Parametrical studies based on the [Formula: see text]-SPH model regarding the breaking and non-breaking waves generated by the flow past a submerged foil are also carried out.

Evolution of Velocity Field and Vortex Structure during Run-Down of Solitary Wave over Very Steep Beach

An experimental results on the spatio-temporal variation of velocity field and vortex structure, generated from the separated boundary layers on the offshore side of the still-water shoreline, during the run-down process of non-breaking solitary waves over a 1:3 sloping beach are presented. Three waves having the incident wave-height to water-depth ratios (H0/h0) of 0.363, 0.263, and 0.171 were generated in a wave flume. Two flow visualization techniques and high-speed particle image velocimetry were employed. The primary topics and new findings are: (1) Mechanism of the incipient flow separation, accompanied by formation of the separated shear layer from the beach surface, is elucidated under the adverse pressure gradient, using the fine data of velocity measurements very close to the sloping boundary. (2) Occurrence of hydraulic jump subsequently followed by development of the tongue-shaped free surface and projecting jet is demonstrated through spatio-temporal variation in the Froude number. It is confirmed by a change in the Froude number from supercritical to subcritical range as the free surface rapidly rises from the onshore to offshore side. (3) A complete evolution of the primary vortex structure (including the core position, vortex size, and velocity distribution passing through the vortex core) is first introduced systematically, together with the illustration of temporal variation in the topological structure. The non-dimensional shoreward distance of the vortex core section decreases with the increase in the non-dimensional time. However, the non-dimensional size height of the primary vortex increases with increasing non-dimensional time. (4) Two universal similarity profiles for both the wall jet flow and the shear layer flow demonstrate independency of the two similarity profiles of the wave-height to water-depth ratio and the beach slope. The similarity profiles indicate the promising collapse of the data from three previous studies for 1:20, 1:10, and 1:5 sloping beaches.

Characteristics of higher-harmonic breaking wave forces and secondary load cycles on a single vertical circular cylinder at different Froude numbers

A 3D numerical two-phase flow model based on solving Unsteady Reynolds-averaged Navier-Stokes (URANS) equations has been used to simulate spilling breaking waves past a single vertical cylinder installed at the edge of a 1:10 slope. The volume of fluid (VOF) method is employed to capture the free surface, and the k−ω Shear-Stress Transport (k−ω SST) turbulence model is used to simulate turbulence effects. Mesh and time-step refinement studies have been conducted by examining the free surface elevations in front of the cylinder and the total horizontal wave forces on the cylinder. These two quantities have also been compared with experimental data from Irschik et al. [21] to validate the present numerical model. The present numerical results are in good agreement with the measured data. The process of breaking waves past the cylinder is described and explained. Moreover, the ‘secondary load cycle’ phenomenon which may lead to higher-harmonic structural response is observed. The characteristics of higher-harmonic wave forces are discussed further at different Froude numbers by changing cylindrical diameters and incident wave properties for both breaking and non-breaking wave cases. For non-breaking wave cases, pronounced secondary load cycles are observed for the Froude number exceeding about 0.4, and the magnitudes of the secondary loads are about 8.1%–12.6% of the peak-to-peak wave force (crest to trough value). However, for breaking wave cases, the critical value of the Froude number to observe the pronounced secondary load cycles reduces to 0.35. The relative magnitude of the secondary load to the peak-to-peak wave force reduces to a range of 2.4%–7.6% due to the high slamming force. The duration of the secondary load cycle is not greater than 0.242T for all the cases, where T is incident wave period. For a vertical cylinder with a larger diameter, the pronounced secondary load cycle occurs with a larger wave steepness.

A Model for the Wind-Driven Current in the Wavy Oceanic Surface Layer: Apparent Friction Velocity Reduction and Roughness Length Enhancement

AbstractA simple analytical model is developed for the current induced by the wind and modified by surface wind waves in the oceanic surface layer, based on a first-order turbulence closure and including the effect of a vortex force representing the Stokes drift of the waves. The shear stress is partitioned between a component due to shear in the current, which is reduced at low turbulent Langmuir number Lat, and a wave-induced component, which decays over a depth proportional to the dominant wavelength λw. The model reproduces the apparent reduction of the friction velocity and enhancement of the roughness length estimated from current profiles, detected in a number of studies. These effects are predicted to intensify as Lat decreases and are entirely attributed to nonbreaking surface waves. The current profile becomes flatter for low Lat owing to a smaller fraction of the total shear stress being supported by the current shear. Comparisons with the comprehensive dataset provided by the laboratory experiments of Cheung and Street show encouraging agreement, with the current speed normalized by the friction velocity decreasing as Lat decreases and λw increases if the model is adjusted to reflect the effects of a full wave spectrum on the intensity and depth of penetration of the wave-induced stress. A version of the model where the shear stress decreases to zero over a depth consistent with the measurements accurately predicts the surface current speed. These results contribute toward developing physically based momentum flux parameterizations for the wave-affected boundary layer in ocean circulation models.

La turbulencia asociada con velocidades orbitales de olas que aún no rompen

Ocean-atmosphere exchange processes are known to decisively determine the sea state, the weather and our planet´s climate. With the ultimate goal of a better understanding of the processes contributing with turbulent kinetic energy into both boundary layers above and below the sea surface, we approach the ocean surface wave phenomenon, and in particular we study the non-breaking waves potential effect. Therefore, measurement of particle velocities were made in a fluid under non-breaking waves, for the purpose of detecting turbulence and its association with the wave steepness. A total of 184 experiments were analyzed, each one with a duration of 3.5 minutes and the presence of approximately monochromatic waves with varying steepness (0.012-0.273). The measurements were carried out in a wave tank with dimensions of 12.26 m × 0.55 m x 0.32 m using an acoustic velocimeter (Vectrino Profiler, Nortek). The u, v and w components of velocity were measured on a 3.5 cm long fluid column in 35 cells of 1 mm in height. During the experiments, the voleocity profile between 1.5 cm and 8.6 cm depth was obtained, referred to the water level in the wave tank (h = 48.6 cm ± 0.5 cm). The waves propagated in deep waters (h/λ > 0.5), where λ is the wave length. Only u(t, z) and w(t, z) components were considered for the analysis since, waves are practical two-dimensional (v(t, z) = 0). Power spectra were calculated in as function of frequency corresponding to u'(t, z) and w'(t, z) turbulent fluctuations, an inertial subrange (isotropic turbulence) was detected in the most of the spectra, for certain depths regardless of the wave steepness. Results from turbulent fluctuations frequency spectra show that eddy size involved in transferring energy to smaller ones, increases with the wave slope.

Parameter estimation of a breaking wave slamming load model using Monte Carlo simulation

For offshore wind turbines (OWTs) located in relatively shallow water, the design is influenced by the occurrence of breaking waves. The strongly nonlinear properties associated with the wave breaking process result in challenges in modelling their impact loads on the structures. The total impact loads are normally calculated as the sum of a slowly varying quasi-static load and an impulsive slamming load. The quasi-static load is normally calculated using Morison’s equation and the slamming load is approximated by the Goda model or the Wienke-Oumeraci model. Given the dynamic properties of OWTs, structural resonances might be excited by the impulsive slamming load. Therefore, there is a clear need to evaluate the response effect excited by the slamming load. In this paper, the response of a vertical pile subjected to a severe breaking wave case is investigated by a combination of data from a large-scale experiment and numerical simulations. The slowly varying quasi-static load obtained in a non-breaking wave packet is modelled using Morison’s equation with the wave kinematics obtained from a fully nonlinear potential flow solver OceanWave3D. The governing parameters used in a slamming load model are estimated using the Monte Carlo method and verified by comparing the experimental data with the numerical simulation results. It is found that the slamming coefficient and the curling factor are close to the values found by the Wienke-Oumeraci model, however the impact duration is significantly larger than the values found by the Goda model and the Wienke-Oumeraci model, which is important for the assessment of the dynamic responses of OWTs.

Lagrangian-like Volume Tracking Paradigm for Mass, Momentum and Energy of Nearshore Tsunamis and Damping Mechanism

There is a gap between model- or theory-based research outputs, which suggest that the runup and amplification of nonbreaking waves generally increase as the sea bottom slopes decrease, and field observations, which indicate that tsunami damage has been rarely reported in places with vast continental shelfs. To resolve this contradiction, we propose a Lagrangian-like volume tracking paradigm to describe the energy, mass, and momentum of travelling nearshore tsunamis and apply the paradigm to analyse the tsunami damping mechanism at typical geophysical scales. The results support the following conclusions: (i) The suggested paradigm is consistent with field observations; continental shelfs with long and mild slopes can effectively diminish tsunami impacts. (ii) Potential energy becomes significant due to the energy transformation process on steeply sloped bathymetries. (iii) On mild-sloped bathymetries, tsunami potential and kinetic energies are conserved until breaking occurs. After breaking, undular bores attenuate tsunami energies effectively. (iv) For extended continental shelf bathymetries, more of the tsunami mass is reflected offshore.

Experimental investigation on non-breaking wave forces and overtopping at the recurved parapets of vertical breakwaters

This paper aims at investigating the intensity of the load applied by non-breaking waves on the recurved parapet wall of vertical breakwaters. For this purpose, 2-D experiments were carried out on a simple vertical breakwater physical model with four different geometries of the parapet (with an “exit angle” of recurve of 0°, 45°, 60° and 90°), under both regular and irregular waves.The tests showed the presence of loads of a quasi-static nature and loads with a typical church roof impulsive behaviour, depending on the (non breaking) wave characteristics and on the shape of the parapet. This impulsive force resulted significantly larger than the quasi-static load of the pure vertical wall case. The maximum measured total force under regular waves was compared with the one obtained by a numerical model based on IHFOAM.For irregular waves with Hs equal to 80% of the freeboard, the maximum measured horizontal load on a 90° recurved wall is equivalent to the force applied by a constant pressure of order 5 ρgHs. It was found, by analysing waves near the leading edge of a wave-group, that history effects are important, i.e., in certain circumstances, a high wave preceded by lower waves generated a larger force than if it had been preceded by waves of the same height. For regular tests, the largest impulsive forces occurred basically at the beginning of the tests, during the initial transient ramp-up period.The irregular wave overtopping was also measured, showing good agreement with the existing predictive formulae and confirming the significant reduction of mean discharge with increasing exit angle and overhanging extension of the parapet.

Numerical study of wave run-up on a fixed and vertical surface-piercing cylinder subjected to regular, non-breaking waves using OpenFOAM

In wave-structure interaction, wave run-up is an important phenomenon that needs to be considered in the design of offshore structures. A thorough understanding of the physics of the nonlinear flow phenomena is necessary for the better insight into the run-up phenomenon. The present work, primarily, is focused on the hydrodynamic simulation of wave run-up and mainly seeks to evaluate the importance of wave scattering types identified by Swanet al. (2005) and lateral progressive edge waves on nonlinear wave amplification around a single fixed cylinder. The analysis is performed numerically using Navier-Stokes solver based on OpenFOAM over a range of wave steepnesses and wavelengths. In order to study wave-cylinder interaction, a numerical wave tank is utilized which includes IHFOAM toolbox for wave generation. For the simulations, the cylinder is assumed vertical and surface piercing with a circular cross-section and the incident wave is considered fifth-order Stokes, regular and non-breaking in deep water. The obtained numerical results that are mainly analyzed in this study include local surface elevation, the scattered wave field around the cylinder and the wave forces. A grid-space and time step refinement study are carried out to evaluate the performance and accuracy of numerical wave tank. The numerical results are compared with experimental data provided by ITTC (OEC), (2013) and a good agreement was obtained which indicates that OpenFOAM is very capable of accurate modeling of nonlinear wave interaction with offshore structures. The harmonic analysis for both short and long wave cases indicates that the scattered wave field around the cylinder involves high harmonics wave run-up. It is confirmed that even the wave-structure interaction caused by the linear incident wave, can lead to weakly nonlinear wave amplification around the cylinder. In the case of long waves cases, the harmonic loading involves at-least third harmonics due to the occurrence of the second wave run-up at the front part of the cylinder.

High frequency resonant response of a monopile in irregular deep water waves

Experiments with a weakly damped monopile, either fixed or free to oscillate, exposed to irregular waves in deep water, obtain the wave-exciting moment and motion response. The nonlinearity and peak wavenumber cover the ranges: $\unicode[STIX]{x1D716}_{P}\sim 0.10{-}0.14$ and $k_{P}r\sim 0.09{-}0.14$ where $\unicode[STIX]{x1D716}_{P}=0.5H_{S}k_{P}$ is an estimate of the spectral wave slope, $H_{S}$ the significant wave height, $k_{P}$ the peak wavenumber and $r$ the cylinder radius. The response and its statistics, expressed in terms of the exceedance probability, are discussed as a function of the resonance frequency, $\unicode[STIX]{x1D714}_{0}$ in the range $\unicode[STIX]{x1D714}_{0}\sim 3{-}5$ times the spectral peak frequency, $\unicode[STIX]{x1D714}_{P}$. For small wave slope, long waves and $\unicode[STIX]{x1D714}_{0}/\unicode[STIX]{x1D714}_{P}=3$, the nonlinear response deviates only very little from its linear counterpart. However, the nonlinearity becomes important for increasing wave slope, wavenumber and resonance frequency ratio. The extreme response events are found in a region where the Keulegan–Carpenter number exceeds $KC>5$, indicating the importance of possible flow separation effects. A similar region is also covered by a Froude number exceeding $Fr>0.4$, pointing to surface gravity wave effects at the scale of the cylinder diameter. Regarding contributions to the higher harmonic forces, different wave load mechanisms are identified, including: (i) wave-exciting inertia forces, a function of the fluid acceleration; (ii) wave slamming due to both non-breaking and breaking wave events; (iii) a secondary load cycle; and (iv) possible drag forces, a function of the fluid velocity. Also, history effects due to the inertia of the moving pile, contribute to the large response events. The ensemble means of the third, fourth and fifth harmonic wave-exciting force components extracted from the irregular wave results are compared to the third harmonic FNV (Faltinsen, Newman and Vinje) theory as well as other available experiments and calculations. The present irregular wave measurements generalize results obtained in deep water regular waves.

On the over-production of turbulence beneath surface waves in Reynolds-averaged Navier–Stokes models

In previous computational fluid dynamics studies of breaking waves, there has been a marked tendency to severely over-estimate turbulence levels, both pre- and post-breaking. This problem is most likely related to the previously described (though not sufficiently well recognized) conditional instability of widely used turbulence models when used to close Reynolds-averaged Navier–Stokes (RANS) equations in regions of nearly potential flow with finite strain, resulting in exponential growth of the turbulent kinetic energy and eddy viscosity. While this problem has been known for nearly 20 years, a suitable and fundamentally sound solution has yet to be developed. In this work it is demonstrated that virtually all commonly used two-equation turbulence closure models are unconditionally, rather than conditionally, unstable in such regions. A new formulation of the $k$–$\unicode[STIX]{x1D714}$ closure is developed which elegantly stabilizes the model in nearly potential flow regions, with modifications remaining passive in sheared flow regions, thus solving this long-standing problem. Computed results involving non-breaking waves demonstrate that the new stabilized closure enables nearly constant form wave propagation over long durations, avoiding the exponential growth of the eddy viscosity and inevitable wave decay exhibited by standard closures. Additional applications on breaking waves demonstrate that the new stabilized model avoids the unphysical generation of pre-breaking turbulence which widely plagues existing closures. The new model is demonstrated to be capable of predicting accurate pre- and post-breaking surface elevations, as well as turbulence and undertow velocity profiles, especially during transition from pre-breaking to the outer surf zone. Results in the inner surf zone are similar to standard closures. Similar methods for formally stabilizing other widely used closure models ($k$–$\unicode[STIX]{x1D714}$ and $k$–$\unicode[STIX]{x1D700}$ variants) are likewise developed, and it is recommended that these be utilized in future RANS simulations of surface waves. (In the above $k$ is the turbulent kinetic energy density, $\unicode[STIX]{x1D714}$ is the specific dissipation rate, and $\unicode[STIX]{x1D700}$ is the dissipation.)