Three-dimensional wave breaking | Nature

Experimental set-up

The experiments were carried out at the FloWave Ocean Energy Research Facility at the University of Edinburgh. The facility has a 2-m-deep, 25-m-diameter, circular wave tank surrounded by 168 flap-type wavemakers. The tank’s circular geometry can create waves travelling in any direction. The wavemakers are operated using a force-feedback control strategy. In this mode of operation, the wavemakers generate and also absorb waves to mitigate the build-up of reflections in the tank.

The surface elevation was measured in the tank using resistive wave gauges. To perform measurements with a high spatial resolution, we developed an 8 × 8 array of wire wave gauges spaced at 0.1 m intervals (array A). The array covers an area of 0.7 m by 0.7 m. Experiments were repeated with the array positioned in six different locations to achieve an effective measurement area spanning from x = −1.4 to 2.8 m and y = 0 to 0.7 m, where (x = 0, y = 0) is the centre of the circular wave tank (Extended Data Fig. 1). We also used a linear array of rigid drop-down wave gauges (array B) to measure the surface elevation over a larger proportion of the tank from x = −8 to 6 m in two locations where y = 0 and y = −0.7 (Extended Data Fig. 1). Measurements from array B were used to confirm that the waves were focused in the region covered by array A, to validate measurements from the newly developed wire wave gauges and, when positioned at y = −0.7 m, to test the symmetry of the waves, which we later assumed when plotting surfaces. The gauges were cleaned and calibrated at the start of each day of testing. Videos of the experiments were recorded using two cameras positioned at the side of the wave tank (Extended Data Fig. 1).

Experimental matrix

To produce breaking waves, we generated steep focused wave groups using linear dispersive focusing. Inputs to the wave tank were defined using linear wave theory. The desired surface elevation

is constructed as a linear summation of wave components with N discrete frequencies that propagate in M discrete directions with frequency ω_n and wavenumber \({k}_{n}=| {{\bf{k}}}_{m,n}| =| [{k}_{n}\cos ({\theta }_{m}),{k}_{n}\sin ({\theta }_{m})]| \) that obey the linear dispersion relation \({\omega }_{n}^{2}=g{k}_{n}\tanh ({k}_{n}h)\) where g is the gravitational acceleration and h is water depth. The direction θ = 0 corresponds to waves that propagate in the positive x direction. z is positive in the upwards direction, with z = 0 corresponding to the still-water level, and t is time. The phases ψ_m,n are defined such that all components are in phase at the desired focus time (t = 0) and position (x = 0, y = 0). Defining the phase in this manner creates a focused wave group, assuming linear dispersive focusing. The duration of each experiment T defines the resolution of the discrete frequencies ω_n = 2πn/T used in equation (1). We set this to 64 s and set the number of discrete directions to M = 144. We left 10 min of settling time between each experiment to allow any unabsorbed reflections in the tank to dissipate, which ensured that each experiment was carried out in as close to quiescent conditions as possible.

We defined the amplitude spectrum of the wave groups we created using the JONSWAP spectrum:

The parametric form of the JONSWAP spectrum in equation (2) generally corresponds to an energy spectrum (and, thus, the amplitude of discrete wave components \({a}_{n}\propto \sqrt{E({f}_{n})}\). Instead, we set the amplitude spectrum to be proportional to the JONSWAP spectrum itself, a_n ∝ E(f_n), as this gives the correct shape of extreme (and, thus, breaking) waves in an underlying random Gaussian sea^38,39. In equation (2), f = ω/2π is frequency, and we set the peak frequency f_p = 0.75 Hz. Here, σ = 0.9 and 0.7 for f < f_p and f > f_p, respectively. We chose a JONSWAP spectrum as this spectral shape represents well typical ocean conditions. We set γ = 1 to give a broad underlying spectrum. A broad spectrum results in a wave group that is well dispersed and less steep at the wavemakers, which minimizes the errors associated with generating linear waves. Moreover, our preliminary experiments found that wave groups with broad underlying spectra exhibited only a single breaking crest, whereas focused wave groups based on narrower spectra were more likely to break several times.

We defined the directional spreading of the wave groups we created using a wrapped normal distribution:

where θ is the angle of propagation, θ₀ is the mean direction and σ_θ is the spreading width. For crossing groups (Δθ ≠ 0), we superimposed two wrapped normal distributions with mean directions θ₀ = ±Δθ/2. When Δθ = 0, we set the mean direction θ₀ = 0. This means that all the wave groups we created had the same mean direction of propagation, which was along y = 0 in the positive x direction. Note that the spreading that we implemented in our experiments was independent of frequency.

Extended Data Table 1 details the different directional distributions we examined. The spectral components a_n corresponding to each experiment were scaled to give the desired input steepness α₀ = a₀k_p at the intended point of linear focus at the centre of the wave tank (x = 0, y = 0), where a₀ is the linearly predicted amplitude at the focus and k_p is the peak wavenumber. We performed several experiments for each directional distribution and varied the input steepness with decreasing increments of Δα₀ reaching a minimum of 0.0125 (Δa₀ ≈ 2 mm) to find the point at which breaking onset occurs. Breaking onset was identified visually. Note that it is also possible to detect breaking using the surface elevation^51,52 and acoustic measurements⁴⁰. Once this threshold was determined, the largest non-breaking wave for each directional distribution was recreated and measured with the high-density gauge array (array A) located in several positions to obtain measurements of surface elevation over the desired area (Extended Data Fig. 1). To understand how directionality affects wave evolution beyond breaking-onset steepness, the same measurement process was also carried out for waves at 112.5%, 125% and 150% of the breaking-onset steepness \({\alpha }_{0}^{\star }\) for selected directional spectra (denoted with a dagger symbol in Extended Data Table 1).

Definitions of the breaking-onset steepness

The quantities used to parameterize wave-breaking onset can have a significant influence on the perceived results. For example, in 2D, the global and local definitions of steepness (see below) can lead to apparently conflicting parameterizations of breaking-onset steepness (as a function of frequency bandwidth)^19,53,54. Additionally, when waves are directionally spread, steepness parameters, which are predominantly 2D, can become ill-defined.

Geometric parameters, such as global steepness or local slope, are generally not considered to be good at distinguishing between breaking and non-breaking waves¹⁸. Despite this, we investigated geometric and spectral measures of steepness for two main reasons. First, steepness-based parameters are the simplest to measure. Second, recent work by ref. ¹⁹ has shown that in 2D, the local slope of waves may function well as a breaking-onset threshold parameter (ref. ¹⁹ also showed that the perceived issues associated with geometric parameters are the result of inconsistent definitions). Although slope may be a useful parameter for indicating the onset of wave breaking, it is a local parameter of individual waves so that it cannot be readily used to predict wave-breaking onset in phase-averaged wave models. Instead, the global steepness can be used to predict the breaking onset within a given sea state. Potentially, it has a broader application for predicting wave breaking in wave forecasting. Thus, we sought to parameterize the breaking onset in terms of both the local slope and the global steepness. By performing experiments with focused wave groups, we could gain a stochastic understanding of how the global steepness may relate to the local wave slope and the wave-breaking onset. This approach is based on the theory of quasi-determinism^38,39, which states that extremes within a sea state exist in the form of wave groups. This relies upon the assumption that the largest waves break. In other words, we are concerned with dominant-wave breaking (at length scales that correspond to the spectral peak).

Global steepness S

The global steepness S provides a linear approximation to the maximum local slope \(\max (\partial {\eta }^{(1)}/\partial x)\) that a wave may have for a given amplitude spectrum (distribution). S is calculated discretely as a sum of N wave components:

The global steepness relates only to the spectrum underlying a given set of waves and does not account for nonlinear wave evolution, the phase coherence (focusing) of individual waves or the directions in which the waves are travelling. Thus, the global steepness can be thought of as a measure of spectral steepness (for amplitude spectra). For focused wave groups, in which the phases of wave components are aligned, the value of S may be realized and exceeded for finite-amplitude nonlinear waves. In conditions where the phases of wave components are not predetermined in such a way or where strong nonlinear focusing occurs, the steepness S may not be related well to the actual slope of the waves in question. Additionally, equation (4) does not take into account the directions in which waves propagate and is a 2D estimate of the slope.

We measured the global steepness of the waves we created in the tank S_M using the frequency spectrum of the surface-elevation time series measured at the centre of the tank. If nonlinear focusing were significant, local changes to the spectral shape would cause the measured value of S_M to differ from the underlying linear value S. The waves produced in the tank were less steep than the input values S₀, and this underproduction by the wavemakers was a function of the overall directional spread (consistent with ref. ⁴¹, which studied less steep, non-breaking, wave groups). As a result, the input values of S₀ were not entirely representative of the waves created in the tank, so instead, we report the measured values S_M. We believe that this decision is justified, owing to the observations we make in Fig. 1d,e, which suggest that for the more directionally spread waves that we created, nonlinear focusing is not significant.

3D ‘global steepness’

We now introduce a measure akin to the global steepness that accounts for the effects of directional spreading. For directionally spread (3D) linear waves, the surface slope has two orthogonal components:

and

The directions in which waves propagate affect how they contribute to the total surface slope. When calculating S in equation (4), which corresponds to the maximum possible 2D linear slope (M = 1 and θ = 0), the phase argument in equation (5) is ignored. A similar approach may be used to calculate a 3D equivalent of S:

For broadbanded spectra, the slope maxima of each component are not necessarily colocated in space. As a result, the 3D global slope S_3D can be quite different from the maximum achievable slope for a linear focused wave group \(\max (| \nabla {\eta }^{(1)}| )\). Thus, we did not use S_3D in ‘Results’. To calculate the linear maximum achievable slope \(\max (| \nabla {\eta }^{(1)}| )({\sigma }_{\theta },\Delta \theta )\), we searched for the time and position at which the absolute slope was a maximum for a focused wave group based on each directional distribution. Extended Data Fig. 2 demonstrates how two focused wave groups with the same global steepness S can have different waveforms and local slopes (Extended Data Fig. 2b,c). The wave group in Extended Data Fig. 2a is unidirectional, whereas the wave group in Extended Data Fig. 2d is an axisymmetric standing wave.

Local steepness

Local measures of steepness, unlike the global steepness, implicitly capture the effects that nonlinear evolution and directionality may have on the waveform of steep waves. As a result, such measures may provide better descriptions of the free surface elevation at breaking onset for any given wave. The local steepness has various forms for which either a locally measured wave amplitude a or a height H (Extended Data Fig. 2) are non-dimensionalized using a chosen length scale, which is most commonly the wavelength λ (ak and kH/2, where k = 2π/λ). The choice of length scale and the manner by which it is calculated (in space or time⁴²) have significant implications for the resulting measure of steepness. In highly spread conditions, even if highly resolved spatial measurements of surface elevation are available, the wavelength and other length scales become ill-defined⁴⁷. Even in 2D, where some of the aforementioned issues are resolved, the local steepness is not always a robust indicator of breaking onset¹⁹. As a result, we have not used local steepness parameters herein.

Local slope ∣∇η∣

The local surface slope is a form of steepness that is well defined and has been shown to be a robust indicator of breaking and non-breaking behaviour in 2D¹⁹. In 3D, the local slope has x and y components. We report the magnitude of the gradient vector \(| \nabla \eta | =\sqrt{{\eta }_{x}^{2}+{\eta }_{y}^{2}}\), where η_x = ∂η/∂x and η_x = ∂η/∂y. The local slope is akin to \(\max (| \nabla {\eta }^{(1)}| )\) but is locally measured not linearly predicted and is, thus, affected by nonlinear wave evolution as well as by directionality.

Measures of spreading for Ω
₀ and Ω
₁

The integral measure of spreading Ω
₀

The integral measure of spreading,

can be used to parameterize the overall spreading of the waves we created, where θ is the direction in which the waves travelled and Ω(θ) is the directional distribution. This parameter is one minus the in-line velocity reduction factor used in ref. ²⁹ and in engineering practice. It is a frequency-independent equivalent of the directional width parameters output by WAVEWATCH III (ref. ⁴³) and ECWAM (ref. ⁴⁴). The parameter Ω₀ provides a measure of the degree to which wave components are standing (as opposed to travelling). In the limits σ_θ → ∞ or Δθ → π, the value of Ω₀ tends to one. It tends to zero for unidirectional waves, σ_θ → 0 and Δθ → 0. The integrand in equation (8) is a function only of θ, as spreading was independent of frequency in our experiments. In the ocean, spreading can be a function of frequency^32,55. Here, we define frequency-independent spreading to reduce the overall complexity of our experiments. For the directional distributions that we define (equation (3)), equation (8) can be expressed as

One limitation of Ω₀ is that it has the same value for two unidirectional counter-propagating wave groups (σ_θ = 0 and Δθ = π) as for an axisymmetric wave group (σ_θ = ∞). These two conditions represent 2D and axisymmetric standing waves, which are known to have different limiting forms^24,45. Figure 3a illustrates how Ω₀ varies as a function of σ_θ and Δθ. The markers show where our experiments are located within this parameter space.

Phase-resolved measure of spreading Ω
₁

As mentioned above, the integral measure of spreading Ω₀ does not fully capture the effects of directionality. We, thus, sought to find an alternative parameter. Directionality affects the shape of wave groups, as demonstrated in Extended Data Fig. 2, which shows the surface elevation (at t = 0) of two wave groups with the same global steepnesses S but different local slopes and steepnesses. Following arguments made in 2D studies, which demonstrated that certain values of local slope ∂η/∂x may trigger breaking^19,54,56, alongside our observations in this paper, which suggest that the focusing of directionally spread wave groups is predominantly linear, we introduced a single-parameter measure of spreading that describes how directionality affects linearly predicted wave slope.

For a given directional spectrum, the degree of spreading affects the maximum surface slope in space and time, which can be predicted linearly as \(\max (| \nabla {\eta }^{(1)}| )\), where η is the free surface surface elevation. Normalizing \(\max (| \nabla {\eta }^{(1)}| )\) by the the 2D global steepness S (\(\max (| \nabla {\eta }^{(1)}| )/S\)) gives a measure of directional spreading, which we will call Ω₁, that reflects how directional spreading affects the potential slope of 3D wave groups (we use the superscript (1) to emphasize that ∇η⁽¹⁾ is predicted linearly). The phase-resolved spreading measure Ω₁ may appear to be a linear measure of slope. However, because we normalize by the global slope (which is a 2D approximation of the slope), \({\varOmega }_{1}=1-\max (| \nabla {\eta }^{(1)}| )/S\) is a measure of directional spreading, which is purely a function of the directional distribution.

Note that, although Ω₀ may not fully describe the effects of directional spreading, it can be calculated simply using operations (integration) on the directional spectrum and is already an output of phase-averaged wave models. Calculating \(\max (| \nabla {\eta }^{(1)}| )\) involves the use of linear wave theory to search for a maximum slope in space and time. In the narrow-banded limit Ω₁ → 1 − S_3D/S, where S_3D is a spectral measure of the slope that ignores phase. Like Ω₀, 1 − S_3D/S is quick to calculate, but its values are quite different to Ω₁ for the broadband spectra in our experiments.

Parametric fitting coefficients

Extended Data Table 2 details the coefficients for the parametric curves fitted to experimentally measured values of local slope ∣∇η∣^⋆ and global steepness \({S}_{{\rm{M}}}^{\star }\) corresponding to breaking onset, which are presented in Fig. 3.

Error quantification

We identified and quantified three main sources of experimental error: wave gauge calibration error, the error associated with the discrete steps by which we varied the input steepness when identifying the breaking onset, and random error, which we estimated from repeated experiments. All three sources of error affected our estimates of the global steepness S_M and the local slope ∣∇η∣. The local slope was further affected by an error that results from estimating the slope from discrete points corresponding to gauges in the high-density gauge array.

Error in the global slope S

First, we will discuss how the above sources of error can affect values of the global steepness at breaking onset (Figs. 2 and 3c,f). The mean calibration error of the wave gauges was ±0.3%, which resulted in an absolute error of Δa₀ = 0.002 m at worst (\(\Delta {S}_{{\rm{M}}}^{\star }=0.008\)). The breaking onset was identified by increasing the input steepness in discrete steps until breaking could be identified visually. We estimated this error as the difference between the threshold values of \({S}_{{\rm{M}}}^{\star }\) (\({\alpha }_{0}={\alpha }_{0}^{\star }\), for the steepest non-breaking waves) and values of S_M calculated for experiments identified as the least steep breaking waves (\({\alpha }_{0}={\alpha }_{0}^{\star }+0.0125\)). The mean value of this error across experiments was \(\Delta {S}_{{\rm{M}}}^{\star }=0.0082\). To quantify the random error, we used the standard deviation of values of \({S}_{{\rm{M}}}^{\star }\) obtained from repeats of the same experiment to quantify the experimental repeatability. The mean value of this standard deviation across our experiments was \(\Delta {S}_{{\rm{M}}}^{\star }=0.011\). If these three errors are treated as independent and combined, the resulting error bars are smaller than the markers used to plot \({S}_{{\rm{M}}}^{\star }\) in Figs. 2 and 3. Therefore, we did not include error bars for measured values of the global steepness at breaking onset \({S}_{{\rm{M}}}^{\star }\) in Figs. 2 and 3c,f.

Error in the post-breaking-onset amplitude a
_M

For the post-breaking-onset behaviour in Fig. 5, the error bars correspond to the standard deviation of the maximum amplitudes measured across repeated experiments. For these experiments, the calibration error was negligible in comparison to the random error obtained from repeated experiments (the error associated with identifying the breaking onset in discrete steps was not applicable to these experiments).

Error in the local slope ∣∇η∣

The high-density gauge array was designed such that the gauges were as closely spaced as possible, while still preventing electric ‘cross-talk’ between the wire gauges so that we could obtain the best possible estimates of local slope. To estimate the local slope ∣∇η∣ (Fig. 3b,e), we performed first-order central differencing, which is associated with a truncation error. To obtain an estimate of this error, we performed a second-order bivariate Taylor-series expansion of η(x, y), from which we obtained an estimate of the error of the gradient vector:

We applied the first-order central differencing twice to obtain estimates of the second derivative and set the error of the magnitude of the local slope to be equal to the magnitude of the error of the gradient vector in equation (10):

For the measured local slope in Fig. 3b,e, the three aforementioned sources of error were negligible in comparison to the truncation error defined in equation (11), and the error bars in Fig. 3b,e correspond to ±Δ∣∇η∣.