Evolution of skewness and kurtosis of weakly nonlinear unidirectional waves over a sloping bottom

We consider the effect of slowly varying depth on the values of skewness and kurtosis of weakly nonlinear irregular waves propagating from deeper to shallower water. It is known that the equilibrium value of kurtosis decreases with decreasing depth for waves propagating on constant depth. Waves propagating over a sloping bottom must continually adjust toward a new equilibrium state. We demonstrate that weakly nonlinear waves may need a considerable horizontal propagation distance in order to adjust to a new shallower environment, therefore the kurtosis can be notably different from the equilibrium value for each corresponding depth both on top of and beyond a bottom slope. A change of depth can provoke a wake-like spatially non-uniform distribution of kurtosis on the lee side of the slope. As an application, we anticipate that the probability of freak waves on or near the edge of the continental shelf may exhibit a rather complicated spatial structure for wave fields entering from deep sea.


Introduction
As waves propagate from deeper to shallower water, their properties are transformed.Classical ray theory predicts slight shortening of wave length, with accompanying change of group velocity, amplitude and steepness, while the frequency is conserved for stationary conditions (Goda, 2010;Mei et al., 2005).The nonlinear stability of wave trains also changes character.In water deeper than the critical depth kh > 1.363, a uniform wave train is modulationally unstable for narrow-band perturbations, while in water shallower than the critical depth kh < 1.363, long-crested waves become modulationally stable (Benjamin, 1967;Whitham, 1974), k being the wavenumber that corresponds to the depth h.
In deep water it is known that third-order nonlinearity can cause focusing of long-crested and narrow-banded waves, and can be responsible for the occurrence of freak waves (Onorato et al., 2001(Onorato et al., , 2005;;Mori and Yasuda, 2002;Janssen, 2003;Mori and Janssen, 2006;Toffoli et al., 2008Toffoli et al., , 2009)).As nonlinear modulation is reduced for long-crested waves on shallower depths, it may be anticipated that there is a corresponding reduction in the probability of freak waves (Mori and Janssen, 2006;Janssen and Onorato, 2007;Toffoli et al., 2009;Janssen, 2009).When the depth is less than the critical depth kh < 1.363, the loss of modulational instability may lead to reduced probability of freak waves (Mori and Janssen, 2006).
Short-crested wave fields on deep water have near-Gaussian statistics despite the nonlinearity of the wave field (Onorato et al., 2002;Socquet-Juglard et al., 2005;Gramstad and Trulsen, 2007;Onorato et al., 2009).In shallower water, however, Toffoli et al. (2009) suggested that when directional components are taken into account, the probability of freak waves is notably increased compared to long-crested waves.In this paper we do not consider the effect of shortcrestedness.
In much shallower water the nonlinear dynamics are dominated by second-order nonlinearity.For this case Pelinovsky and Sergeeva (2006) studied the probability of freak waves with the Korteweg-de Vries equation, and found that the skewness and kurtosis was controlled by the Ursell number.In this paper we limit consideration to depths similar to or deeper than the critical value; third-order rather than secondorder nonlinearity is important for nonlinear evolution.We do however account for second-order static nonlinear contributions in the reconstruction of the wave field.
A number of studies have been dedicated to nonlinear propagation of wave groups over a slope.Djordjević and Redekopp (1978) discussed a solution for an envelope-hole soliton and the fission of an envelope soliton moving over variable depth.Benilov et al. (2005) described how a wave packet propagating from deep to shallow water can be decreased in amplitude, and if the depth reaches the critial value, the packet can disperse.Fructus and Grue (2007) Published by Copernicus Publications on behalf of the European Geosciences Union.
simulated the shoaling of a solitary wave group on a slope in shallow water, and they found that waves steepen with reducing depth.Grimshaw and Annenkov (2011) considered how a wave packet evolves close to the critical depth kh = 1.363.Recent experiments considered wave amplifications and phase modulation over a shoal (Jarry et al., 2011).
Recently the statistics of wavefields over a slope, relevant for freak wave occurrence, have also been studied.Janssen and Herbers (2009a,b) studied the kurtosis of waves with directional spreading over a submerged shoal in deep water.They found that strongly non-Gaussian behavior can be produced due to the concomitant effects of focusing and nonlinearity.Furthermore, the maximum of kurtosis is found downwave of the position where the waves are steepest near the top of the shoal.In these cases, the deeper kh is around 20 and the shallower kh is 0.22.Sergeeva et al. (2011) investigated the transformation of a random wave field in variable shallow water depth through a Korteweg-de Vries framework.They found that kurtosis increases as the water depth decreases and the kurtosis can achieve a local maximum close to the shallower side.In this case, the deeper kh is 0.44 and the shallower kh is 0.3 and the depth gradient on the slope ∂h/∂x = 0.03.Experimental evidence presented by Trulsen et al. (2011) shows that there can be a local maximum of kurtosis and skewness close to the shallower side of the slope.In these experiments, the maximum deeper kh was 1.6, the minimum shallower kh was 0.54 and the depth gradient on the slope was ∂h/∂x = 0.05.
In the present study we limit attention to waves on finite depth (kh ≥ 1.2) on a mild slope.We employ one of the simplest nonlinear models capable of describing this situation, namely the nonlinear Schrödinger equation with variable coefficients and a shoaling term for slowly varying depth, first derived by Djordjević and Redekopp (1978) and subsequently refined by Peregrine (1983) and Iusim and Stiassnie (1985) and enhanced with the effects of a nonuniform current by Turpin et al. (1983).They assumed small wave steepness = ka 1, mild slope ∂h/∂x = O( 2 ) and finite depth (kh) −1 = O(1).More sophisticated models are also available, e.g.Liu and Dingemans (1989) relaxed the requirement of slowly varying depth to ∂h/∂x = O( ), while Xiao and Lo (2004) extended the analysis to higher order effects of variable depth beyond shoaling with the assumption ∂h/∂x = O( 4/3 ).More exact nonlinear models have also been employed in recent studies (Fructus and Grue, 2007;Janssen and Herbers, 2009a,b).Here we shall suffice with the simplest model that includes leading order nonlinearity, dispersion and shoaling.
It is known that the coefficient of the cubic nonlinear term vanishes for kh = 1.363.This fact has prompted several workers to develop higher-order nonlinear terms in the evolution equation applicable to both flat bottom (e.g.Johnson, 1977;Kakutani and Michihiro, 1983;Slunyaev, 2005) and to sloping bottom (Grimshaw and Annenkov, 2011).They suggest a rescaling for kh ≈ 1.363 in order to reestablish a leading order balance between nonlinearity and dispersion.The rescaling implies longer spatial and temporal scales, or equivalently, narrower bandwidth.We take the opposite approach, in the sense that we let the incoming wave fields have non-vanishing bandwidths.Indeed, in an analysis of wave records from the North Sea, a likely region for application of the present results, Trulsen and Dysthe (1997) observed that wave fields of overall steepness O( ) could have a larger bandwidth of O( 1/2 ), thus linear dispersion is expected to dominate over nonlinearity also away from the critical depth kh = 1.363.
Our numerical results reveal that the kurtosis is reduced for smaller depths.We also find that the nonlinear dynamical response of the wave train may occur over a longer domain than the region of depth change, such that the wave train may need a long propagation distance after the slope in order to reach the new equilibrium values for both kurtosis and skewness.Our main result is that the characteristic relaxation distance can be large compared both to the local wavelength and the scale of depth change.As a consequence, we anticipate that waves entering a continental shelf from deeper water can have a spatially non-uniform distribution of freak waves, notably different from that expected from equilibrium statistics on the given depth.

Nonlinear Schrödinger equation for uneven bottom
We assume irrotational flow of an inviscid and incompressible fluid with a free surface.A coordinate system is oriented with horizontal axis x along the quiescent water level toward shallower water, the z-axis vertically upward, and with the bottom located at z = −h(x).The velocity potential is (x,z,t) and the free surface is located at z = η(x,t).The governing equations are and The characteristic frequency ω c of the incident waves and the acceleration of gravity g are used for normalization.A characteristic wavenumber for the deeper water can be estimated by the linear dispersion relation for infinite water depth k c = ω 2 c /g.The surface displacement is characterized by the amplitude a, and we define a small ordering parameter by means of the steepness of the waves = ka 1 which is assumed to be small.We also assume that the depth is finite but not too small, (k c h) is slowly varying, ∂h/∂x = O( 2 ), and that the modulation of the waves in space and time is slow on scales x and t.
In the following, all variables and parameters are scaled and made dimensionless according to these assumptions.Throughout the following, all reference to space and time, x and t, imply the slow modulation space and temporal variables mentioned above.
After normalization, the characteristic frequency is 1, and the characteristic wavenumber is the solution of 1 = ktanh(kh).We introduce the following perturbation series for the surface displacement η and velocity potential where θ = x k(ξ )dξ − t, and c.c. denotes complex conjugate.
The equations of Djordjević and Redekopp (1978) are expressed in terms of the velocity potential, but here we express the equations in terms of the first harmonic amplitude B of the surface displacement.By sorting orders and harmonics, we get the following result for the first harmonic at the third order where c g is the group velocity and σ = tanh(kh).We also get the following result for the zeroth harmonic at the third order In these equations the four quantities h, k, σ and c g are all functions of x.
We take the following general solution for the latter equation, and substitute into the former equation and obtain the nonlinear Schrödinger equation with shoaling term for uneven bottom The four coefficients are In our case we let f (x) = 0 because the term with f (x) only affects the phase of the wave (see Djordjević and Redekopp, 1978).From the plot of the first three coefficients as a function of dimensionless depth in Fig. 1, it is seen that the coefficient of the shoaling term is much smaller than the coefficients of the nonlinear term and the dispersion term.Figure 1 also shows that in shallow water the nonlinear coefficient varies rapidly for small depth.From the derivation, we also get Finally, the reconstruction of surface elevation η is achieved by invoking Eq. ( 5).

Model setup
In the following we present results for one bathymetry.The region x ≤ 0 has a constant depth h = h 1 , while the region x ≥ L has a constant depth h = h 2 with h 2 < h 1 .In the region x ∈ (0,L), there is a slope defined by h = h 1 + (h 2 − h 1 ) x L .The bathymetry is shown in Fig. 2.
The wave field is assumed to be periodic in time with period T .The temporal Fourier transform is where t n = n t and j = j ω with j ranging from 1 − N/2 to N/2.We let the length of each time series be T = 1000 which corresponds to approximately 159 periods.We let N = 320, thus the time step is t = T N = 3.125 and ω = 2π T ≈ 0.0063.Unidirectional incoming waves with an initially Gaussian spectrum have been studied.The Fourier amplitudes at the starting point x = x 0 are given by where the phases ϒ j are statistically independent and uniformly distributed on the interval [0,2π).We let = 0.1.The Benjamin-Feir index is defined by Janssen (2003) in terms of frequency bandwidth as where we compute the bandwidth δ using the spectral halfwidth at the half maximum, δ = σ ω √ 2ln2.Simulations with σ ω = 0.05, 0.085, 0.17 which correspond to the initial values for BFI =1.7, 1, 0.5 have been carried out.To solve the spatial evolution of the wave field according to Eq. ( 10), we have used Matlab's ode45, which is a fourth order Runge-Kutta solver with variable step.Due to the application of periodic boundary conditions in time, FFTroutines are used to calculate the temporal derivatives.In the simulation, the accuracy of the solution is controlled by the pre-chosen relative and absolute error tolerances in the solver.These have been chosen so that the error is less than 10 −3 .

Results
Monte-Carlo simulations have been performed to investigate the evolution of skewness and kurtosis over the sloping bottom.The statistical parameters are computed as functions of location.An example of the statistical convergence of kurtosis with respect to ensemble size is shown in Fig. 3.An ensemble size of 3000 seems to be sufficient for our simulations, and has been chosen in the following.
In order to ensure it is a realistic wave field that first arrives at the slope, we insist that the value of kurtosis should have reached a stable value before the wave field reaches the slope.Figure 4 shows a typical evolution of kurtosis on a flat bottom after initialization at x 0 = 0.
In the following we initialize at some x 0 < 0 sufficiently far from the beginning of the bottom slope so that a stabilized wave field has been established before reaching the slope.
Figure 5 is a contour plot of stabilized kurtosis with respect to different initial BFI and 1/(kh) after imposing the initial condition on flat bottom.Cubic interpolation was used to plot the contours.When the depth kh is greater than 1.363, high initial BFI and large depth lead to high stabilized kurtosis.However, when kh = 1.363, nonlinear effects disappear and the kurtosis is equal to 3.Moreover, the result agrees well with previous publications that stabilized kurtosis is larger in deep water and smaller in shallow water (Janssen, 2003;Janssen and Onorato, 2007).
We have considered the following cases of finite water depths h 2 =1, 1.1955, 2, 3, 4 corresponding to k 2 h 2 =1.2, 1.363, 2.065, 3.015, 4.003.The deep water depth has been set to h 1 = 10, corresponding to k 1 h 1 = 10.In all cases we let the bottom slope be located between x = 0 and x = L = 200.The kurtosis and skewness of the surface elevation with three different initial BFI are shown in Figs.6-11.It is shown that when the initial BFI is not less than 1, there is a local minimum of kurtosis near the shallower side of the slope.Furthermore, for the same initial BFI, the natural response length for k 2 h 2 = 1.2 is always the longest among all cases.When k 2 h 2 > 1.363, the skewness has the same behavior as the kurtosis.Nevertheless, when k 2 h 2 ≤ 1.363, the skewness achieves a minimum over the slope and then relaxes toward the value at the equilibrium state on the shallower side of the slope.Furthermore, the skewness for the case k 2 h 2 = 1.2 has the "sharpest" minimum at the end of the slope among all cases.This could be provoked by the sudden change in the bottom slope.For all cases, when the shallower dimensonless water depth is below the threshold k 2 h 2 < 1.363, the dynamic reponse happens over a much longer distance than the bottom slope itself and it takes quite a long relaxation length for both the kurtosis and skewness to become stabilized in the shallow domain.
Figure 12 shows kurtosis variation with respect to k 2 h 2 after the slope for different initial BFI.As we see from the figure, with diminishing depth, the kurtosis decreases for all initial BFI.When k 2 h 2 is larger than 1.363, cases of larger initial BFI will have larger stable kurtosis.However, when k 2 h 2 is less than 1.363, cases of smaller initial BFI get larger stable kurtosis after the slope.Comparing Fig. 12 with Fig. 5, it can be found with the condition that the depth of a flat bottom is the same as the depth of the shallower side x  for a sloping bottom, and with the same initial BFI for incoming waves, the stabilized kurtosis after a sloping bottom is smaller than for a flat bottom when k 2 h 2 > 1.363, larger when k 2 h 2 < 1.363 and the same when k 2 h 2 = 1.363.We have also carried out simulations with several other bathymetries in addition to the one shown here.Those additional simulations indicate that when the depth gradient is large on the shallow side of the slope, there can be a long natural response length.However, when the depth gradient is small on the shallow side of the slope, kurtosis and skewness can stabilize almost at the same location as the change of depth.Moreover, when the shallower depth is less than the critical depth kh < 1.363, the wave field always needs a certain relaxation length to become stabilized in the shallower domain.

Conclusions
Through Monte-Carlo simulations using a nonlinear Schrödinger equation for variable depth, we have investigated some effects of nonuniform bathymetry on wave statistics.We have found evidence that a local change in bathymetry may provoke non-equilibrium statistics in a region that may extend far beyond the local depth nonuniformity.We have found evidence that when a wave field travels over a bottom slope into shallower water, a wake-like structure may be anticipated on the shallower side for the skewness and the kurtosis.As an application, we anticipate that the probability of freak waves on or near the edge of the continental shelf may exhibit a rather complicated spatial structure for wave fields entering from deep sea.

Fig. 3 .
Fig. 3. Typical convergence of the ensemble-averaged kurtosis with respect to the number of runs in the ensemble.

Fig. 12 .
Fig. 12. Kurtosis variation with respect to k 2 h 2 after the slope for different initial BFI.