Jie Zhang

Michel Benoit
email: benoit@irphe.univ-mrs.fr

Yuxiang Ma

Equilibration process of out-of-equilibrium sea-states induced by strong depth variation: Evolution of coastal wave spectrum and representative parameters

Keywords:

Recent studies showed both experimental and numerical evidence that the occurrence probability of freak waves could be significantly enhanced as results of non-equilibrium dynamics induced by strong depth variations. The sea-state is characterized by strong non-Gaussian behavior in a short spatial extent after the depth transition, covering a few wavelengths. In this work, we investigate the complete equilibration process of an out-of-equilibrium sea-state via high-fidelity numerical simulations. In the simulations, the region after the depth transition is set as long as around one hundred wavelengths, such that the spectral adaptation develops and terminates eventually. The results are analyzed with spectral, cross-spectral and statistical approaches. It is shown that there are two stages with different spatial scales in the equilibration process. In the short scale, the sea-state is characterized by significant changes in wave statistics, freak wave occurrence probability is intensified. In the long scale, the wave spectrum undergoes strong modulation, the spectral peak disintegrate into a relative broad band, and low-frequency waves are enhanced as well. We show evidence that the spectral changes in the long scale are due to interactions of free components. The wave nonlinearity is shown to be positively correlated to the magnitude of the dynamical responses, but irrelevant to the length of the spatial scales in the equilibration process. In the established shallow-water equilibrium, the freak wave occurrence probability becomes less than Gaussian expectation and the waves are asymmetric in the vertical direction and symmetric in the horizontal.

Introduction

In oceanography, the term "freak wave" (or "rogue wave") denotes a wave whose crest-to-trough height takes a very large value in comparison to the background waves. One common criterion used to identify a freak wave is > 2 , where is the significant wave height of the sea-state (see e.g., [START_REF] Draper | Freak' ocean waves[END_REF][START_REF] Haver | Evidences of the existence of freak waves[END_REF]. Such extreme waves have attracted considerable attention of the scientific community, since they pose severe risks to the safety of maritime structures, ship navigation and human lives. In recent decades, more and more freak wave events were observed and reported (see e.g., Nikolkina and Didenkulova, 2011b;[START_REF] O'brien | Extreme wave events in Ireland: 2012-2016[END_REF][START_REF] Didenkulova | Freak waves in 2011-2018[END_REF][START_REF] Häfner | Real-world rogue wave probabilities[END_REF], implying that freak waves may occur with startling frequency in the real world and potentially manifest under a broad range of conditions (Nikolkina and Didenkulova, 2011a). Combining the statistics in Nikolkina and Didenkulova (2011b) and [START_REF] Didenkulova | Freak waves in 2011-2018[END_REF], plotted in Fig. 1, it is seen that from 2006 to 2018, over 80% of the freak waves occurred either in the shallow water zone (with water depth lower than 50 m following their definition) or at the coast, causing tragic loss of lives. Moreover, the ships and infrastructures damage induced by freak waves (not listed here) took place mostly in the shallow and coastal regions too. Thus, it is of great importance to better evaluate the risk of freak waves, especially in coastal areas.

The topic of freak waves has been studied for several decades, with a number of hypotheses put forward to explain the occurrence of abnormal large waves from both linear and nonlinear perspectives (see e.g., [START_REF] Dysthe | Oceanic rogue waves[END_REF][START_REF] Kharif | Rogue Waves in the Ocean[END_REF][START_REF] Adcock | The physics of anomalous ('rogue') ocean waves[END_REF][START_REF] Dudley | Rogue waves and analogies in optics and oceanography[END_REF]. However, the scientific community has not yet come to a consensus on the dominant mechanism of freak waves occurring in nature. The modulation instability as a result of four wave quasi-resonant interactions [START_REF] Benjamin | The disintegration of wave trains on deep water Part 1. theory[END_REF] has gained much attention and is proposed as the prototype of freak waves [START_REF] Onorato | Extreme waves, modulational instability and second order theory: wave flume experiments on irregular waves[END_REF][START_REF] Akhmediev | Waves that appear from nowhere and disappear without a trace[END_REF][START_REF] Wang | Experimental study of breathers and rogue waves generated by random waves over non-uniform bathymetry[END_REF], because it allows extremely large waves to manifest spontaneously in a local scale. However, some researchers argue that the wave modulation instability can hardly survive under complex conditions in real seas [START_REF] Fedele | Real world ocean rogue waves explained without the modulational instability[END_REF]. For example, the modulation instability mechanism fails to explain the freak waves occurring in the coastal areas with relative water depth ℎ ( denotes the wavenumber and ℎ the local water depth) lower than 1.363, since the modulation instability ceases to exist. Recent studies revisited the well-known Draupner wave and attributed its formation to the particular sea conditions, including the significant wave height, the spectral shape, and the directional condition [START_REF] Cavaleri | The Draupner wave: A fresh look and the emerging view[END_REF]. The Draupner wave has been reproduced in a wave basin with good quality by separating the sea-state into two crossing systems [START_REF] Mcallister | Laboratory recreation of the Draupner wave and the role of breaking in crossing seas[END_REF]. This work indicates that the limiting crest height due to wave breaking could be significantly elevated in crossing sea-states.

In the coastal areas, significant decrease of water depth is another popular mechanism of freak waves. [START_REF] Trulsen | Laboratory evidence of freak waves provoked by non-uniform bathymetry[END_REF] showed experimental evidence that, as unidirectional irregular waves passed over a steep slope, strong non-Gaussian behavior of the sea-state occurred within a short range after the end of the slope, and the occurrence probability of extreme waves was significantly enhanced. The underlying physics is attributed to the non-equilibrium dynamics (NED) [START_REF] Trulsen | Rogue waves in the ocean, the role of modulational instability, and abrupt changes of environmental conditions that can provoke non equilibrium wave dynamics[END_REF]: due to depth variations, the sea-state leaves the quasi-equilibrium status it had in the deeper region and adapts to a new shallow-water equilibrium.

The investigation of the depth variation induced NED effects has been recently enriched by a number of experimental works. For waves propagating over submerged bars, [START_REF] Ma | Experimental study of statistics of random waves propagating over a bar[END_REF][START_REF] Ma | Variations of statistics for random waves propagating over a bar[END_REF] showed the spatial evolution of statistical parameters of the sea-state, and the relation between skewness and kurtosis of the free surface elevation (i.e., normalized third and fourth order statistical moments). [START_REF] Chen | Effect of bottom slope on the nonlinear triad interactions in shallow water[END_REF] considered different up-slope gradients of a submerged bar, and studied the strong nonlinear triad interactions occurring during the shoaling process using wavelet bispectral analysis. Recent experimental work reported in [START_REF] Trulsen | Extreme wave statistics of long-crested irregular waves over a shoal[END_REF] (hereafter referred to as TRJR20 for brevity) showed that the condition ℎ < 1.3 on the relative water depth is required after the depth transition, such that the NED enhances the occurrence probability of freak waves conspicuously. For waves propagating over shoals, [START_REF] Kashima | Estimation of freak wave occurrence from deep to shallow water regions[END_REF] showed the maximum wave heights as functions of skewness and kurtosis. [START_REF] Kashima | Aftereffect of high-order nonlinearity on extreme wave occurrence from deep to intermediate water[END_REF] attributed the occurrence of freak waves to the aftereffects of modulation instability.

In the large-scale experiments reported in [START_REF] Zhang | Statistics of extreme waves in coastal waters: large scale experiments and advanced numerical simulations[END_REF], transient second-order harmonics appeared atop of the slope, and the bound harmonics were shown to play an important role. Experiments of waves passing over abrupt depth changes were also conducted. [START_REF] Bolles | Anomalous wave statistics induced by abrupt depth change[END_REF] showed that, in the out-of-equilibrium area, the distribution of free surface elevation conforms closely to a Gamma distribution. Based on this distribution, the skewness, kurtosis and other statistical properties can be estimated.

The NED effects were also studied with numerical simulations. [START_REF] Janssen | Nonlinear wave statistics in a focal zone[END_REF] was the pioneering work pointing out that the NED could elevate the probability of freak waves basing on Monte-Carlo simulations with a spectral wave model. [START_REF] Sergeeva | Nonlinear random wave field in shallow water: variable Korteweg-de Vries framework[END_REF] then studied the NED effects with the Korteweg-de Vries model, showing the increase of skewness and kurtosis with the decrease of water depth. [START_REF] Zeng | Evolution of skewness and kurtosis of weakly nonlinear unidirectional waves over a sloping bottom[END_REF] adopted Monte-Carlo simulations with the nonlinear Schrödinger equation, and found a complex spatial variation of wave statistics over the bottom step. [START_REF] Gramstad | Freak waves in weakly nonlinear unidirectional wave trains over a sloping bottom in shallow water[END_REF] tested the NED effects resulting from different shapes of bottom profile using a Boussinesq model with improved dispersion. [START_REF] Viotti | Extreme waves induced by strong depth transitions: fully nonlinear results[END_REF] firstly investigated the problem using fully nonlinear potential flow (FNPF) simulations (basing on the free surface Euler equations), and attributed the non-Gaussian behavior of the sea-state to the spectral settling effect. The effects of directionality turn out to be significant as well: [START_REF] Ma | A numerical investigation on nonlinear transformation of obliquely incident random waves on plane sloping bottoms[END_REF] showed large values of the incident angle of a wave train weaken triad-wave interactions (with a fully nonlinear Boussinesq-type model), and [START_REF] Ducrozet | Influence of varying bathymetry in rogue wave occurrence within unidirectional and directional sea-states[END_REF] showed a broader angular spreading of wave energy suppresses the occurrence of freak waves using simulations with a model based on the high-order spectral (HOS) method. [START_REF] Zheng | Fully nonlinear simulations of unidirectional extreme waves provoked by strong depth transitions: The effect of slope[END_REF] studied the experiments of [START_REF] Trulsen | Laboratory evidence of freak waves provoked by non-uniform bathymetry[END_REF] with a FNPF model solved by the boundary element method, and considered the roles played by different harmonic components. Their results implied that the secondorder harmonics dominated the increase of freak wave probability. [START_REF] Zhang | Wave-bottom interaction and extreme wave statistics due to shoaling and de-shoaling of irregular long-crested wave trains over steep seabed changes[END_REF] (hereafter referred to as ZB21) further investigated the experiments of TRJR20 with a different FNPF model, called Whispers3D. In their work, the importance of the third-order harmonics on the freak wave occurrence probability was revealed via the harmonic separation technique [START_REF] Fitzgerald | Phase manipulation and the harmonic components of ringing forces on a surface-piercing column[END_REF].

From the theoretical perspective, Didenkulova and Pelinovsky (2016) showed the dispersion focusing conditions for freak wave formation in strongly inhomogeneous media in linear regime, while [START_REF] Onorato | Twenty years of progresses in oceanic rogue waves: the role played by weakly nonlinear models[END_REF] showed, in weakly nonlinear regime, that the skewness and kurtosis are related to the change of spectral width. Very recently, Li et al. (2021b) demonstrated with second-order asymptotic theory that, after abrupt depth transitions, the peaks of statistical parameters arise from the interaction of free components and bound components in a wave packet. Such an explanation is in line with the conclusion drawn in ZB21 for irregular wave trains. Later, their second-order theory was also verified by numerical simulations (Li et al., 2021c) and experiments (Li et al., 2021a).

In summary, the NED could result in strong non-Gaussian behaviors of the sea-state, including: (i) local variations of the statistical parameters of the free surface elevation and the kinematics (orbital velocities, accelerations, etc.) beneath the free surface, (ii) spectral modulations due to rapid wave energy exchanges among harmonics, and (iii) a heavy tail of the wave height exceedance probability distribution. Previous studies have shown that the NED exists in a longer spatial extent than that of the depth variation, the delay of dynamics being called "latency" in [START_REF] Zheng | Fully nonlinear simulations of unidirectional extreme waves provoked by strong depth transitions: The effect of slope[END_REF]. The spatial desynchronization between the depth variation and the change of wave dynamics poses risks of freak waves not only in the sloping bottom area, but also in a certain region after it. The knowledge on the spatial scale where the NED affects the sea-state is therefore of practical interest. However, most of the existing works focus on the area shortly after the depth variation (in a few wavelengths), and rarely discuss the aftereffects of the NED in the wave evolution in a larger scale. This work extends the analysis of the step-bottom case discussed in ZB21 to a much longer spatial extent, from ( ) to (10 2 ) in the shallow water area, with being the wavelength of the peak period. Here, the equilibration process of the sea-state is investigated indepth, and the wave evolution is characterized by combining spectral and statistical analyses to study the transformation of the wave spectrum and statistical characteristics of the wave field.

The remainder of this paper is organized as follows: in section 2, we introduce the considered cases, seabed configurations and incident wave conditions. Section 3 shows the essence of the adopted FNPF model used for wave modelling as well as the selected numerical parameters, while section 4 presents the methods adopted for interpreting the results. In section 5, the complete process of the sea-state equilibration is analyzed, and the effects of wave nonlinearity on the wave evolution are discussed. The main findings and conclusions are summarized in section 6.

Bathymetric configuration and wave conditions

In the experimental work of TRJR20, the non-equilibrium phenomenon has been studied with unidirectional irregular waves propagating over a submerged bar. The bar is of steep slopes on two sides (±1∕3.81), both provoking strong seastate adaptation processes. To determine the condition of the occurrence of NED, TRJR20 conducted a series of tests with a similar level of wave nonlinearity but different relative water depths. They showed that the condition ℎ < 1.3 over the shallower region was required in their tests for the manifestation of strong non-Gaussian behavior of the sea-state. The Run 3 in TRJR20 is of particular interest, as it contains measurements of the horizontal flow velocity below the still water level (SWL) on one hand, and as the kurtosis of the free surface elevation achieves a maximum value as high as 4.2 (indicating a heavy tail of the wave height distribution) on the other hand. In ZB21, deterministic numerical simulations were performed and validated by comparing with the measurements of TRJR20, showing the high-fidelity of the code. Then, additional tests with a step-bottom profile (not considered in TRJR20) were simulated to isolate the NED effects induced by the shoaling and de-shoaling processes.

Here, we pursue the investigation of the step-bottom case introduced in ZB21, by extending the shallower region far after the step, over a distance of order (10 2 ), to study the complete stabilization process of the out-of-equilibrium sea-state induced by the depth transition.

The considered setup and the numerical wave tank are shown in Fig. 2 (note the origin of -axis is placed at the start of the shallower region). Similar to ZB21, the bottom profile consists of a deeper flat zone with uniform water depth ℎ 1 = 0.53 m and a shallower zone with ℎ 2 = 0.11 m, connected by a steep (1/3.81) plane slope. The peak wave period is fixed for all cases considered here:

= 1.1 s. The corresponding wavelengths are ,1 = 1.80 m in the deeper region and ,2 = 1.07 m in the shallower region. ,2 is used hereafter to normalize the distance in the shallower region. The length of the shallower region is set either to 103.6 m in case 1 or to 43.6 m in the other cases, which correspond to roughly 96 ,2 and 40 ,2 , respectively.

In ZB21, the Run3 of TRJR20 was simulated in a deterministic way, i.e. we used the record of free surface elevation at the first probe located 1.1 m before the depth transition in the experiment of Run 3 to force the numerical model. Furthermore, the simulations lasted for the same duration as in the experiment (about 4,900 waves with period ). Here, a different approach is adopted: the incident irregular wave trains imposed at 0 = -2.7 m are generated by linear superposition (as described in Appendix A) from a JONSWAP spectrum, governed by a significant wave height , a peak frequency = 1∕ and a peak enhancement factor :

( ) = 2 (2 ) 4 1 5 exp - 5 4 4 exp -( -) 2 2 ( ) 2 , ( 1
)
where denotes the gravitational acceleration, is the adjustment factor for , and is the spectral asymmetry parameter ( = 0.07 for ≤ and = 0.09 for > ). The significant wave height is defined here as = 4 √ 0 , with being the ℎ -order moment of a given spectrum

( ): = ∫ ∞ 0 ( )d .
(2)

The JONSWAP parameters of the reference case 1 are set as = 0.030 m, = 1.1 s and = 3.3. Six additional cases are introduced to investigate the effects of the incident sea-state nonlinearity, by varying from 0.010 m to 0.035 m. The configurations of all tested cases are summarized in Table 1, where three non-dimensional parameters are computed for both the deeper and shallower zones, including relative water depth ≡ ℎ, wave steepness ≡ with ≡ √ 2 0 , and Ursell number ≡ ∕ 3 . a For all cases, = 1.1 s, = 3.3. Parameters with subscript 1 refer to the deeper region with depth ℎ 1 = 0.53 m, and those with subscript 2 refer to the shallower region with depth ℎ 2 = 0.11 m.

Numerical model and parameters

Wave propagation model

Wave propagation is modelled here using the FNPF theory, assuming (i) the fluid is inviscid and incompressible, (ii) the flow is irrotational, and (iii) the surface tension is negligible. Consider a two-dimensional Cartesian coordinate system ( , ) with the origin set on the SWL and -axis pointing upward, the computational domain is bounded by a moving free surface = ( , ) and a fixed impermeable bottom profile = -ℎ( ).

The equations governing wave motion are

+ = 0 in -ℎ ≤ ≤ , (3)
+ - = 0 on = , (4) + 1 2 2 + 2 + = 0 on = , ( 5
) ℎ + = 0 on = -ℎ, (6)
where ( , , ) represents the velocity potential, and subscripts denote partial derivatives. Following [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], the free surface boundary conditions ( 4) and ( 5) can be expressed as functions of the free surface variables and ̃ ≡ ( , = , )

= -̃ + ̃ 1 + 2 , ( 7
) ̃ = -- 1 2 ̃ 2 + 1 2 ̃ 2 1 + 2 , ( 8
)
where ̃ ( , ) ≡ ( , = , ) is the vertical component of the velocity on the free surface. ̃ is needed to integrate Eqs. (7-8) in time, determining it from the known quantities ( , ̃ ) is the so-called Dirichlet-to-Neumann (DtN) problem. The DtN problem is solved in different models such as Boussinesq-type models [START_REF] Madsen | A Boussinesq-type method for fully nonlinear waves interacting with a rapidly varying bathymetry[END_REF][START_REF] Bingham | Velocity potential formulations of highly accurate Boussinesq-type models[END_REF], HOS models [START_REF] Dommermuth | The initialization of nonlinear waves using an adjustment scheme[END_REF][START_REF] Gouin | Development and validation of a non-linear spectral model for water waves over variable depth[END_REF] and the Hamiltonian Coupled-Mode model (Belibassakis and Athanassoulis, 2011; [START_REF] Papoutsellis | Implementation of a fully nonlinear Hamiltonian coupled-mode theory, and application to solitary wave problems over bathymetry[END_REF] with different methodologies.

In the recent work of [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF], the DtN problem is solved by using a Chebyshev-tau approach initially introduced in [START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF]. With this method, the FNPF problem is solved without perturbation expansion (thus no need to evaluate high-order derivatives), and no additional assumption made on the wave nonlinearity, dispersion, nor on the bathymetry profile. This model, called Whispers3D, is therefore considered as fully nonlinear and dispersive. It has been presented and validated in previous works [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF][START_REF] Benoit | Analysis of the linear version of a highly dispersive potential water wave model using a spectral approach in the vertical[END_REF][START_REF] Zhang | Statistics of extreme waves in coastal waters: large scale experiments and advanced numerical simulations[END_REF][START_REF] Zhang | Wave-bottom interaction and extreme wave statistics due to shoaling and de-shoaling of irregular long-crested wave trains over steep seabed changes[END_REF], to which interested readers are referred for a detailed description. Only the main steps are outlined below.

First, the physical domain with varying boundaries ∈

[-ℎ, ] is mapped into the computation domain ∈ [-1, 1]
with fixed boundaries, via the change of vertical coordinate:

( , , ) = 2 + ℎ - ℎ + . ( 9
)
Then, the governing equations (3-6) are reformulated in the transformed ( , , ) domain, and we denote the potential with ( , ( , , ), ) ≡ ( , , ). Secondly, the potential is approximated using a set of orthogonal Chebyshev polynomials of the first kind, denoted as ( = 0, 1, ..., ):

( , , ) ≈ ( , , ) = ∑ =0 ( , ) ( ), (10)
where the coefficients ( , ), = 0, 1, ..., , are now the main unknowns of the problem. Thirdly, the Chebyshevtau method is applied to the Laplace equation ( 3), providing -1 equations. The problem is closed with the bottom boundary condition ( 6) and a Dirichlet free surface boundary condition ( , = 1, ) = ̃ . Thus, a linear system with + 1 equations for + 1 unknowns is built. Lastly, with the coefficients determined from the linear system, the vertical velocity at the free surface is computed as

̃ ( , = 1, ) = 2 ℎ + ∑ =1 2 , ( 11
)
so that Eqs. (7-8) can be integrated in time.

It should be mentioned that, although effects like viscous bulk dissipation and wave breaking are readily included in Whispers3D by modifying the free surface boundary conditions Eqs. (7-8) [START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF], they are not considered in the present work. Such a choice allows simplifying the discussion and to focus on the effects resulting from NED. The overall dissipation would decrease the level of sea-state nonlinearity as waves propagate, and the wave breaking effects would possibly limit the height of the highest waves. It is anticipated that such effects would introduce some differences, but not overturn our conclusions on the NED effects.

Numerical parameters and settings

In Whispers3D, Eqs. (7-8) are marched in time using an explicit strong-stability preserving third-order Runge-Kutta scheme [START_REF] Gottlieb | On high order strong stability preserving Runge-Kutta and multi step time discretizations[END_REF], with a constant time-step Δ . Spatial derivatives are approximated using high-order finite difference schemes on a regular grid with spatial resolution Δ using centered stencils of nodes. The accuracy and efficiency of the model can be balanced by properly choosing the numerical parameters: , , Δ and Δ . Following the successful simulation of the Run 3 of TRJR20 reported in ZB21, the same set of numerical parameters is adopted here, namely = 7, = 5, Δ = 0.01 m, and Δ = 0.01 s. These time and spatial steps correspond to roughly 110 points per peak wave period and peak wave length ,2 . The corresponding Courant-Friedrichs-Lewy numbers are CFL 1 = ,1 Δ ∕( Δ ) = 1.64 in the deeper region and CFL 2 = 0.97 in the shallower one.

In the numerical model, waves are generated and damped using the relaxation zone technique [START_REF] Bingham | A Fourier-Boussinesq method for nonlinear water waves[END_REF], see the setup of the relaxation zones in Fig. 2. Methods used to discretize the target spectrum and construct the free surface elevation and potential of the incident wave field by linear superposition are classic, they are briefly re-called in Appendix A.1. Absorbing long waves in the relaxation damping zone revealed challenging during the simulations. Therefore, a specific choice of the damping zone parameters (in particular its length ) was made, see Appendix A.2.

Numerical wave gauges are evenly distributed along the flume, recording the free surface elevation with sampling frequency = 1∕Δ = 100 Hz. For ∈ [-2.6, 3] m, the spacing between wave gauges is 0.2 m, and for ∈ [3, ] m, it is 0.5 m.

To obtain a sufficient number of waves for subsequent spectral and statistical analyses and to evaluate the uncer-tainty of the simulated results, each case of interest is simulated over a duration of 500 s, equivalent to 455 waves (with period ) and repeated 10 times with different random phases in the imposed time series. Eventually, for each case, the size of the data-set is about 4, 500 waves.

Methods of sea-states analysis

Statistical and spectral parameters

The NED is induced by rapid changes of water depth, and the sea-state adaptation process is characterized by strong non-Gaussian behavior. To investigate the equilibration process of an out-of-equilibrium sea-state, we start by computing statistical parameters and applying a spectral (Fourier) analysis technique.

In statistical analysis, we consider high-order moments of the free surface elevation. The mean value ⟨ ⟩ (with ⟨⋅⟩ denoting an averaging operator) is removed from the computed signal . Then, the standard deviation = √ 0 = √ ⟨ 2 ⟩ can fully describe the statistical properties of a Gaussian sea-state. For characterizing non-Gaussian characteristics of an out-of-equilibrium sea-state, the normalized thirdand fourth-order moment, namely skewness 3 and kurtosis 4 , are needed:

3 = ⟨ 3 ⟩ 3 , 4 = ⟨ 4 ⟩ 4 . ( 12
)
3 indicates the overall vertical asymmetry of the signal (i.e. with respect to the horizontal SWL). In a Gaussian sea-state, 3 = 0 is expected, indicating that the signal is symmetric with respect to the SWL. When nonlinearity plays a role, 3 takes non-zero values. Positive skewness means waves have sharp high crests and shallow flat troughs in general, and vice versa for negative values. 4 is related to the number of large values, and thus the extremes in the signal. In a Gaussian sea-state, 4 is exactly equal to 3, and larger values of 4 indicate an increased occurrence of extreme events compared to the Gaussian (linear) expectation.

The asymmetry parameter, , is another parameter that characterizes the non-Gaussianity of a sea-state. It is defined on the basis of bispectrum, it can also be computed as the skewness of the Hilbert transform of , denoted ( ) (Elgar and Guza, 1985a):

= ⟨ ( ) -⟨( )⟩ 3 ⟩ ⟨ ( ) -⟨( )⟩ 2 ⟩ 3∕2 , ( 13
)
offers a measure of the overall horizontal asymmetry of the signal (i.e. with respect to a vertical axis): a positive means that waves are tilting forward as a whole, while negative values indicate waves tilting backward.

In spectral analysis, the variance density spectrum is evaluated with the Welch method. The time series is divided into segments with 8, 192 points (approximately 82 s ≈ 75 ) and 50 % overlapping rate between every two consecutive segments. Such choices result in a relatively fine spectral resolution Δ ≈ ∕75 ≈ 0.012 Hz.

For a given spectrum ( ), the spectral width is characterized by the peakedness factor (see [START_REF] Goda | Random seas and design of maritime structures, Third Edition[END_REF]:

= 2 2 0 ∫ 2 ( )d . (14)
Here, following e.g. [START_REF] Waseda | Evolution of a random directional wave and freak wave occurrence[END_REF], we consider and plot 1∕ , since it increases as the spectrum broadens. The peak frequency is estimated as

̂ = ∫ 4 ( )d ∫ 4 ( )d . ( 15
)
In Eqs. ( 14) and ( 15), = 0.5 and = 0.5 are set to discard the low-frequency (LF) components ( ∈ [0, 0.5 ]) in the spectrum while evaluating the spectral parameters.

Cross-spectral analysis

It is known that Fourier analysis assumes a stationary and ergodic sea-state and no coupling among components. In the current scenario, the NED effects clearly violate such assumptions. It brings some inconvenience while interpreting the Fourier spectrum, since the free and bound harmonics cannot be distinguished without further information. To achieve a better understanding of the physical processes taking place over the submerged step, researchers have adopted different approaches. The bi-spectral analysis technique [START_REF] Hasselmann | Bispectrum of ocean waves[END_REF] has been used to detect triad nonlinear interactions (Elgar and Guza, 1985a;[START_REF] Ma | Experimental study of statistics of random waves propagating over a bar[END_REF][START_REF] Chen | Effect of bottom slope on the nonlinear triad interactions in shallow water[END_REF][START_REF] Zhang | Statistics of extreme waves in coastal waters: large scale experiments and advanced numerical simulations[END_REF][START_REF] Zhang | Wave-bottom interaction and extreme wave statistics due to shoaling and de-shoaling of irregular long-crested wave trains over steep seabed changes[END_REF]. However, such a method does not provide information about the ratio between free and bound components.

The harmonic separation technique [START_REF] Baldock | A laboratory study of nonlinear surface waves on water[END_REF][START_REF] Fitzgerald | Phase manipulation and the harmonic components of ringing forces on a surface-piercing column[END_REF] is a powerful tool to isolate the time series that correspond to harmonics at different orders with few repetitions of the tests using specific phase shifts in the incident wave field. It has been applied in analyzing the bound harmonics that appear atop of the bottom slope, and provided valuable insights [START_REF] Zheng | Fully nonlinear simulations of unidirectional extreme waves provoked by strong depth transitions: The effect of slope[END_REF][START_REF] Zhang | Wave-bottom interaction and extreme wave statistics due to shoaling and de-shoaling of irregular long-crested wave trains over steep seabed changes[END_REF]. However, as this method is built on the basis of Stokes expansions, its applicability to uneven bottom cases remains to be further investigated.

The cross-spectral analysis has been used to study the dispersive properties of both deep water waves and shoaling waves (see [START_REF] Mitsuyasu | On the dispersion relation of random gravity waves. Part 2. an experiment[END_REF][START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF]Elgar and Guza, 1985b). Recently, [START_REF] Martins | Dispersive characteristics of non-linear waves propagating and breaking over a mildly sloping laboratory beach[END_REF] studied strongly nonlinear waves propagating over a mild slope with a surf zone using cross-spectral analysis. This method is applied in our study and briefly described hereafter.

We start by computing the cross-spectrum 1 , 2 ( ) from two time series recorded at nearby locations 1 and 2 , using similar Welch method and settings as presented in the previous sub-section. The spectra of phase lag 1 , 2 ( ) and coherence coh

1 , 2 ( ) can then be obtained as

1 , 2 ( ) = arctan ⎡ ⎢ ⎢ ⎢ ⎣ ℑ 1 , 2 ( ) ℜ 1 , 2 ( ) ⎤ ⎥ ⎥ ⎥ ⎦ , ( 16
) coh 1 , 2 ( ) = 1 , 2 ( ) 1 , 2 ( ) * 1 , 1 ( ) 2 , 2 ( ) 1∕2 , ( 17
)
where ℑ(⋅) and ℜ(⋅) denote the imaginary and real parts of a complex number, respectively. Here, the distance between

1 and 2 separating two consecutive wave probes is either 0.47 ,2 , or 0.93 ,2 depending on the position.

The phase lag 1 , 2 ( ) is bounded between -∕2 and ∕2. By unwrapping the phase lag 1 , 2 ( ), the real phase change Θ

1 , 2 ( ) can be recovered for each wave component. The unwrapped phase lag allows determining the spectrum of the dominant wavenumber ( ):

( ) = Θ 1 , 2 ( ) 2 -1 , ( 18
)
which measures the combined effects of free and bound components at ( 1 + 2 )∕2. For a free wave with frequency , its wavenumber follows the dispersion relationship, and the relative water depth is

tanh ( ) = ℎ, ( 19
)
where = 2 ∕ . If a wave with frequency is a n ℎorder bound super-harmonic, its wavenumber follows the dispersion relationship of the primary harmonic ( ∕ ), and the relative water depth follows

tanh ( ) = ℎ∕ 2 . ( 20
)
For a wave with frequency being a n ℎ -order bound subharmonic, then its relative water depth is

tanh ( ) = 2 ℎ. ( 21
)
ℎ can be viewed as an averaged measure of the relative water depth. When the sea-state is composed of free components only, ℎ is expected to be superimposed with , whereas when bound harmonics exist, ℎ deviates from . Hereafter, we discuss deviations of ℎ from toward 2 , 2 and 3 , and assume higher order effects being negligible in the considered cases.

Analysis of the sea-state equilibration process

In this section, the various parameters, spectra and crossspectra shown in the figures represent the mean of 10 realizations of each case (with different sets of phases). The 95% confidence intervals (CI) are calculated for the parameters ̂ , 1∕ , , 3 , and 4 .

Long scale sea-state evolution

We first show the long scale sea-state evolution by analyzing case 1. The shallower flat region is 103.6 m (≈ 96 ,2 ) long in this case, such that the full process of the sea-state equilibration process can be observed.

Fig. 3 shows the spatial evolution of wave spectrum, which is normalized as ( ) ∕ 0 for the convenience of comparison among cases with different levels of energy (recalling = 1∕1.1 Hz denotes the peak frequency of the target JON-SWAP spectrum). ̂ ( ) evaluated over the range [ , ] is superimposed to track the evolution of the main spectral peak. It is seen that, shortly after the abrupt depth change, a particular pattern appears near 2 . This spatial pattern is due to the interactions of the second-order free and bound components, which has also been observed in TRJR20 and explained in ZB21 on the basis of considerations from Massel (1983). Such a trend remains visible, though progressively attenuated over a distance of about 10 ,2 after the depth variation.

As waves propagate further over the shallower region, the spectrum undergoes more pronounced changes. The spectral peak "disintegrates" into a relatively broad band near , whereas the spectral evolution near 2 becomes negligible. Meanwhile, the LF waves are considerably enhanced over the shallower region. The disintegration of the spectral peak develops as waves propagate relatively far from the depth variation. Such a large scale spectral evolution was also observed in the experiments of [START_REF] Zhang | Statistics of extreme waves in coastal waters: large scale experiments and advanced numerical simulations[END_REF], but has not yet been thoroughly analyzed.

In Figs. 4 and5, more detailed spectra at eight specific positions are shown in both logarithmic and linear scales to better highlight the spectral changes. Fig. 4 illustrates the short-scale evolution from the wavemaker ( = -2.7 m) to a position close to the end of the slope ( = 1.6 m). Figs. 5 shows the long-scale evolution from = 10 m till = 100 m, where the spectral peak disintegration process is included. In each panel of two figures, the JONSWAP spectrum is superimposed for comparison. In Fig. 4(a), the target JONSWAP spectrum at the wavemaker is superimposed for comparison. For probes over the shallower region, in Figs. 4(b-d) and 5, the incident JONSWAP spectrum is amended to the shallower water depth ℎ 2 as ( ) ℎ ( ), where ℎ is a shoaling factor from the linear wave theory, computed as

ℎ ( ) = (ℎ 1 , )∕ (ℎ 2 , ), (22)
with (ℎ, ) = d ∕d denoting the group velocity following the linear dispersion relationship. The JONSWAP spectrum corrected with ℎ thus represents the linear expectation of the spectral evolution from the deeper region to the shallower. In Fig. 4(a), we see only long waves with frequencies around 0.1 are enhanced very slightly. The good agreement between the simulated and the target spectrum indicates the validity of the linear wave-making approach and relaxation absorption zone setup. In Figs. 4(b-d), the frequency range in < 1.5 is still in good agreement with JONSWAP spectrum with shoaling coefficient, whereas the second-and third-harmonics manifest and evolve as waves propagate in the area closely after the end of the slope. As can be seen from both Fig. 3 and Fig. 4(b) however, right at the end of the sloping area at = 0, the second order peak around 2 is hardly visible, and the wave spectrum is very close to the incident one, apart from a higher level of LF energy. The weak enhancement of super-harmonics at the end of the slope is due to the fact that this slope is steep and over a short distance, its length (1.6 m) representing 89% of the incident peak wavelength ,1 . Therefore, the dynamical responses develop slower compared to the depth transition.

As waves propagate further over the shallower region, it is observed in Fig. 5 that the energy of high-order harmonics returns to its initial level, the main spectral peak pro- In all panels, the black line represents the spectrum from the numerical simulations and the red line represents the target JONSWAP spectrum amended to water depth ℎ 2 by applying a shoaling coefficient from the linear wave theory.

gressively decreases in magnitude and evolves into a broader band. The spectral shape in the frequency range < 1.5 significantly differs from linear expectation, the modified JON-SWAP spectral shape. Besides, the spatial evolution of the high-frequency (HF) range ( > 1.5 ) is very limited, the HF spectral tail decays algebraically with the same rate as in the modified JONSWAP spectrum, close to -5 . As a whole, Figs. 3-5 clearly show how the wave spectrum "reacts" to the abrupt depth change, first in the HF range with energy transfers towards super-harmonics, and then at a longer scale around the spectral peak and LF range.

The spectral parameters are computed to quantify the spectral changes over the shallower flat region. In Fig. 6(a), the averaged spectral peak frequency ̂ is normalized by of the incident sea-state. We verify ̂ is equal to at the wavemaker, and remains nearly unchanged shortly after the slope. More marked changes occur between = 5 ,2 and = 35 ,2 , the peak frequency shifts first downwards then upwards, and eventually tends to a constant mean level slightly higher than the incident . Note that the most evident frequency downshift takes place at ≈ 10 ,2 , which can be identified in the spectrum shown in Fig. 5(a). The spectral width parameter 1∕ is shown in Fig. 6(b). A small spectral broadening takes place right after the slope, which is related to the enhancement of the higher order harmonics. Then, the spectrum significantly broadens in the range between = 5 ,2 and = 35 ,2 with the evolution of ̂ , this is due to the spectral peak disintegration. Eventually, 1∕ is almost doubled from 0.28 at the incidence to 0.51 at = 35 ,2 .

In Fig. 6(c), the significant wave heights computed for the LF range ∈ [0, 0.5 ], the short-wave range ∈ [0.5 , 0.5 ] and the total range ∈ [0, 0.5 ] of the wave spectrum are shown. They are normalized by the corresponding significant wave heights , of the incident JONSWAP spectra, listed in Table 1. It is noticed that, over the shallower region, of the LF components increases at the expense of of the short-wave components and stabilizes after = 35 ,2 . As the dissipation due to friction or breaking events are excluded, of the total spectral range remains more or less unchanged in the computation domain.

The statistical moments are shown in Fig. 7. In Fig. 7(a), the skewness 3 increases significantly in a short region after the slope and achieves its local maximum (slightly larger than 1) at ≈ 0.75 ,2 . It indicates that the overall wave profile evolves from nearly symmetric to skewed shape with sharp crests and flat troughs. Previous works have explained this local increase as the results of the development of highorder harmonics, and the interactions of bound and free components. After this significant enhancement, 3 gradually decreases to a mean constant level at ≈ 35 ,2 . After that position, 3 ≈ 0.57 indicates that waves remain asymmetric, due to the nonlinear shallow water effects. The evolution of the asymmetry parameter in Fig. 7(b) shows that after entering the shallower region, the waves first tilt backward then forward. It is noticed that, = 0 at = 0.75 ,2 , where 3 achieves its local maximum value. It indicates that the wave profile changes independently in the horizontal and vertical directions, and that the waves are symmetric in the horizontal direction while the most asymmetric profile is achieved in the vertical.

The evolution trend of kurtosis in Fig. 7(c) is similar to that of 3 : a considerable enhancement of 4 takes place in the shallower region close to the end of the slope, indicating that the occurrence probability of extreme waves is elevated. A maximum value of 4.6 is achieved, at the same position ( ≈ 0.75 ,2 ) as for 3 ; then a mild decrease develops in a relatively long scale and 4 eventually converges to a constant level slightly lower than 3, the risk of freak waves is then lower than linear expectation. All parameters shown in Figs. 6 and 7 converge to constant levels roughly at the same position ( ≈ 35 ,2 ), after which the sea-state achieves a quasi-equilibrium. It indicates that the adaptation of a sea-state from a deeper to a shallower regime finishes around = 35 ,2 in case 1. During the equilibration process of the out-of-equilibrium sea-state, two different spatial scales can be identified over the shallower region: a short scale ( ) where the statistical parameters vary significantly, and the probability of freak waves is strongly enhanced; a long scale (10 ∼ 10 2 ) where the spectral shape modification is non-trivial, and the freak wave occurrence probability gets suppressed. The short-scale parameter variations are mainly due to the effects of higher order harmonics, whereas the long-scale variations are due to the spectral disintegration that occurs around the spectral peak.

Role of nonlinearity on the sea-state equilibration

As seen in the previous section, the NED affects the seastate in a relatively long spatial extent, influencing the statistics of the wave field in the area close to the end of the slope, and the spectral shape within a long shallow water area. The non-equilibrium dynamical responses of a sea-state result from both external factors such as bottom gradient, depth change, and internal factors likes relative water depth, seastate nonlinearity, and spectral width. Among the internal factors, the effective relative water depth was proven important, since it dominates the occurrence of NED. Here, the effects of sea-state nonlinearity are investigated. To this end, cases 2-7 are set with increasing significant wave heights (see Table 1). As shown above, the shallow water equilibrium is achieved around = 35 ,2 , thus the shallower region length is reduced to 40 ,2 in cases 2-7 to limit the computational burden.

The evolutions of the normalized spectra of cases 2-7 are shown in Fig. 8 panels (a-f) correspondingly. In case 2 shown in Fig. 8(a), the sea-state is of the mildest nonlinearity. The beating pattern in the frequency range ∈ [1.5 , 2.5 ] is rather clear, and lasts far after the slope. Over the shallower region, few energy is transferred to LF components. The primary components in ∈ [0.5 , 1.5 ] are more or less unchanged. This is because low sea-state nonlinearity after the slope results in limited wave-wave interactions, and thus limited energy redistribution.

When the incident sea-state nonlinearity is increased, the spectral evolution becomes more complex. In the short scale ( ), the spectral changes in all six cases are similar, apart from more energetic second-and third-order harmonics in the beating patterns. Then, in the long scale (10 ∼ 40 ), the effects of sea-state nonlinearity on the spectral evolution become important. As , of the incident sea-state increases in cases 2-7, the disintegration of the spectral peak becomes more and more pronounced, resulting in increasingly significant broadening of the spectrum in the frequency band [0.5 , 1.5 ]. In the most nonlinear situation of case 7, Fig. 8(f) shows that the spectral peak almost flattens in the shallow water equilibrium. The energy of HF harmonics with > 1.5 increases with incident sea-state nonlinearity, with more marked maximum values in the beating area of the 2 peak. In the meantime, on the contrary to what one could have expected, it is observed that the beat- ing pattern sustains shorter and shorter in space. In the cases with relatively strong nonlinearity, the beating pattern is interrupted by the strong nonlinear wave-wave interaction after the waves propagate over a distance of a few wavelengths in the shallower region. Compared to case 2 in Fig. 8(a), stronger incident sea-state nonlinearity results in higher LF wave energy over the shallower region, confirming that the LF waves are generated as results of wave-wave interaction, with limited effect of reflection from the damping relaxation zone.

Fig. 9 shows the evolutions of the spectral parameters in cases 2-7. It is observed in Fig. 9(a) that the peak frequency ̂ ≈ and shows negligible spatial modulation in case 2 with the lowest nonlinearity. The short scale variation of ̂ shows barely no relevance to , of the incident sea-state, since the curves of ̂ ∕ corresponding to cases 2-7 are almost superimposed for < 3 ,2 . Distinct spatial modulations of ̂ take place in the long scale, a frequency downshift first takes place, and followed by an evident upshift before the shallow water equilibrium is achieved. Such a trend of peak frequency evolution gets more marked with increasing , .

In Fig. 9(b), it is shown that the spectral width increases in both spatial scales, and the spectral broadening becomes more pronounced with larger incident , . In case 2, the spectral width is nearly unchanged along the domain, whereas in the cases with higher nonlinearity, the value of 1∕ shows a small raise in the shorter scale, then a significant increase in the longer one. The significant wave height evolution is shown in Fig. 9(c). As the dissipation is not included in the simulations, the variation of total defined in ∈ [0, 0.5 ] is negligible. Only corresponding to LF waves in < 0.5 (dash lines) and short-waves > 0.5 (solid lines) are plotted in Fig. 9(c). The level of LF wave energy clearly increases for higher incident nonlinearity. Such a trend was anticipated, since stronger wave-wave interaction results in more energy transfer from the primary modes to LF modes.

To further illustrate the relation between the spectral parameters and the incident nonlinearity of incident sea-state, Fig. 10 shows the spectral parameters in the equilibrium state (evaluated by taking the arithmetic mean for > 35 ,2 ) as functions of incident ,1 given in Table 1. The Ursell number ,1 is chosen to highlight that the NED of the sea-state is dominated by the combined effects of relative water depth and the nonlinearity. In the present study, the relative water depth of the incident sea-state is fixed, such that the change of ,1 is purely due to the variation of 1 i.e., the nonlinearity. From Fig. 10(a), it is seen that in the cases with relatively mild nonlinearity ( ,1 < 4×10 -3 ), the peak frequency ̂ in the equilibrium is lower than , , and as the nonlinearity of the incident sea-state increases, the ̂ in the steady seastate gradually shifts upward. In Fig. 10 with ,1 in a quasi-linear manner. The variations of for both LF and short wave components as functions of ,1 are shown in Fig. 10(c). The value of for LF components increases, whereas the values of of short waves declines as ,1 increases. The quasi-linear trends are also clear for both two curves.

The statistical moments of the free surface elevation for cases 2-7 are shown in Fig. 11. Evidently, the incident seastate nonlinearity significantly affects the magnitude of the statistical parameters, 3 , and 4 , in their short scale evolutions. In the long scale, where the statistical parameters undergo mild decreases, the effects of nonlinearity on different parameters are dissimilar. In Fig. 11(a), 3 achieves increasingly higher levels in the shallow water steady state with the growth of , . This is again explained by nonlinear effects on shallow water waves, which result in asymmetric wave profiles with sharp high crests and flat shallow troughs in general. However, for the values of and 4 , the relevance of the established levels in the equilibrium zone and the incident sea-state nonlinearity seems less marked. In all cases shown in Fig. 11(b), tends to a constant value, slightly larger than 0. It indicates that the wave symmetry in the horizontal direction is restored as the sea-state adaptation finishes, and that the horizontal asymmetry of waves mainly results from the NED. In Fig. 11(c), 4 of cases 2-7 all converge to a mean level lower than 3, indicating the occurrence probability of freak waves returns close to Gaussian expectation. Such a decrease of freak wave occurrence probability is almost independent of what the level of wave nonlinearity the sea-state had before the slope, and how much 4 was enhanced in upstream.

It is illustrated in Fig. 7 and 11 that significant variations of statistical parameters take place shortly after the depth variation before achieving their steady levels in the long scale. Therefore, not only the steady levels of 3 , and 4 established in the shallow water equilibrium, but also their extreme values achieved in the short scale evolution are shown in Fig. 12. We note that the short scale extremes (maximum and minimum values) of all three parameters increase almost linearly as the incident sea-state nonlinearity augments. Besides, the steady level of 3 shown in Fig. 12(a) is of a linear dependence with ,1 as well. As observed in Fig. 11, the steady levels of and 4 seem to be independent on ,1 . As general remarks of Figs. 9-12, the incident nonlinearity dominates the importance of non-equilibrium dynamical response of the sea-state, the spectral properties in the long scale as well as the statistical parameters in the short scale increase with ,1 in a quasi-linear manner. Whereas, the steady levels of and 4 achieved in the long scale, the positions corresponding to local extremes of statistical parameters, and length required for the sea-state to restore shallow water equilibration are seen to be independent on the incident nonlinearity. Note that such a remark for nonlinearity is obtained with other effects like breaking events and bottom friction excluded. In real sea-states, such a linear trend may not be as straightforward as it is in the idealized cases, but we conjecture that such effects will not overturn the essence of our conclusions. The cross-spectral analysis results are given in Fig. 13 and 14. The spectrum of relative water depth of the mean wavenumber and the spectrum of coherence are shown. The relative water depth for linear harmonics and bound harmonics 2 , 2 and 3 are superimposed for comparison. In the coherence spectrum, we adopt the same thresh-old coh 2 1 , 2 = 0.25 (horizontal dash line) as in [START_REF] Martins | Dispersive characteristics of non-linear waves propagating and breaking over a mildly sloping laboratory beach[END_REF] for analyzing the obtained results.

In Fig. 13, the spectra of ℎ and coh 2 1 , 2 for four specific positions in the short scale of sea-state equilibration is shown. Three representative cases are considered: the quasilinear case 2, the weakly nonlinear case 4 and the highly nonlinear case 7. In the short scale, the spectral modulation takes place mainly over the HF range. We therefore limit our discussion up to third-order harmonics < 3.5 . The fourth-and higher-order harmonics are omitted, since they are of low level of energy, thus with minor contribution to wave evolution. It is noticed that in Fig. 13(a.1), the LF components in the spectrum of ℎ is superimposed with secondorder sub-harmonics 2 . This is due to the fact that the spectral components with < 0.45 are cut-off from the imposed target spectrum, thus waves in this frequency range are bound components resulting from nonlinear wave-wave interactions. Waves with > 0.45 are well described by linear dispersion, since they are imposed as free components. As waves pass the slope, show that the enhancement of second-and third-order bound harmonics is positively related to wave nonlinearity, since the curve of ℎ deviates further from for larger , . The strong deviation of the mean relative water depth ℎ from to 2 and 3 indicates the increased portion of bound waves for a given frequency. It is also noticed that the LF waves for < 0.45 are gradually released as free components in the area from -2.7 m to 1.6 m. The coherence of the two adjacent probes remains high (coh 2 1 , 2 > 0.25) for < 3.5 in Figs. 13(a.2-d.2), except in case 4, low coherence results in a crude estimate of ℎ for 3 < < 4 .

The long scale evolution of the cross-spectral parameters is shown in Fig. 14. Given the spectral modulation mainly occurs around the spectral peak, we discuss the results for < 2.5 . Higher harmonics with > 2.5 are of low level of energy as illustrated in Fig. 8 for ∈ [5 ,2 , 35 ,2 ]. Furthermore, the estimate of ℎ is crude, since the coherence is low for < 2.5 . As a general remark, all three cases with different levels of nonlinearity basically follow the linear dispersion relationship, including the LF components. Only ℎ of case 7 with the strongest incident nonlinearity slightly deviates from for ∈ [2 , 2.5 ]. Such an observation indicates the spectral evolution in the long scale results from the interactions of free components, which makes an important distinction between short scale and long scale wave evolution.

Conclusion

We performed fully nonlinear numerical simulations of unidirectional irregular waves propagating over a submerged step, aiming at the investigation of evolution characteristics of an out-of-equilibrium sea-state. The tested cases are variations of existing experiments reported in TRJR20, and the adopted numerical code has been validated thoroughly in ZB21. The depth change forces the incident sea-state to leave the equilibrium it had in the deeper region and evolve to a new shallow-water equilibrium in a certain distance. While reestablishing the new equilibrium, the sea-state is significantly affected by NED effects. The present work is distinguished from others since the full equilibration process is discussed for the first time, by setting a sufficiently extended shallower region after the depth variation.

It is found that the equilibration process of an out-ofequilibrium sea-state has two spatial scales after the depth transition area. The short scale refers to the area a few wavelengths O( ,2 ) closely after the slope. In this scale, the NED effects result in rapid energy transfer to high-order harmonics and evident change of wave statistics, in particular a significant increase of freak wave occurrence probability. Despite affecting the spectrum in the high-frequency range, the NED effects on the representative spectral parameters are limited. Besides, the cross-spectral analysis shows higher incident nonlinearity results in increased portion of bound super-harmonics in the short scale.

The longer scale covers an extent of dozens of wavelengths O(10 ∼ 10 2

,2 ) after the end of the slope. In the longer scale, the NED leads to remarkable spectral changes in the primary and LF wave components for < 1.5 . The "disintegration" of the spectral peak takes place in the long scale sea-state evolution: the wave energy is redistributed in the vicinity of the peak frequency, resulting in considerable spectral broadening. In the meantime, the LF waves are also enhanced. The mechanism behind the spectral change in the long scale evidently differs from that in the short scale, as the cross-spectral analysis shows that most of the waves evolving in the long scale behave as free harmonics. Although strong spectral modulations, the changes of wave statistics are very limited, highlighting the importance of higher harmonics on the enhancement of freak wave probability.

After the two-stage equilibration process, all parameters converge to constant mean levels correspondingly. In the established shallow-water steady state, the wave spectrum is characterized by a broader band with less pronounced peak, the general wave profile is asymmetric in the vertical direction with sharp high crests and flat shallow troughs. In addition, kurtosis being lower than 3 in the steady state implies that freak wave probability reduces to a level lower than the Gaussian expectation, i.e. a lighter tail of the wave height distribution is expected.

The incident wave nonlinearity was found to dominate the magnitude of the dynamical response of sea-state, including the extremes of statistical parameters in the short scale, and the change of spectral parameters in the long scale.

But it has minor effects on the length of the two scales in the equilibration process, and the positions where the statistical and spectral parameters achieve their maximum values in two scales are not changed for different incident nonlinearity. It should also be noticed that, changing the incident wave nonlinearity will not affect the levels of 4 (lower than 3) and (close to 0) in the steady state ( > 35 ,2 here). From the practical perspective, the findings of the present study would be of interest in the design of coastal structures (site selection) and in the optimization of ship routes. also the LF components. This length corresponds to three wavelengths of a LF wave with = 0.15 Hz or two wavelengths of a LF wave with = 0.1 Hz. Although a bit costly, such a choice is capable of effectively reducing most of the long-wave reflection.

The performance of the damping zone was assessed by calculating the averaged reflection coefficient ̄ , defined as

̄ = ∫ ( ) ( )d ∫ ∕2 0 ( )d , (A.6)
where ( ) is the spectral reflection coefficient, and control the spectral range for the evaluation of ̄ . With the present choice of relaxation zone function, the total wave reflection coefficient, computed for = 0 and = 0.5 , ̄ is around 4% before the slope in all simulations. Furthermore, the reflection of LF waves ( ̄ evaluated for = 0 and = 0.4 ) is very low, varying from 0.01% to 0.05% depending on the nonlinearity level of the case.

Figure 1 :

 1 Figure 1: Statistics of the freak waves reported in mass media over the period 2006 -2018, based on the works of Nikolkina and Didenkulova (2011b) and Didenkulova and Pelinovsky (2020). The number and percentage of freak wave events are shown in (a), and the loss of lives in (b).

Figure 2 :

 2 Figure 2: Definition sketch of the bathymetry adopted in the simulations, with denoting the length of the shallower region, and gray zones representing the relaxation zones.

Figure 3 :

 3 Figure 3: Case 1: spatial evolution of normalized wave spectrum, log 10 ( ( ) ∕ 0 ), plotted in the colormap with a logarithmic scale. The close-to-horizontal dash white line represents the evolution of the peak frequency ̂ , and the vertical ones outline the extent of the bottom slope.

Figure 4 :Figure 5 :

 45 Figure4: Case 1: normalized wave spectra ( ) ∕ 0 at different locations close to the depth transition. In all panels, the black line represents the spectrum from the numerical simulations. In panel (a) the blue line is the target JONSWAP spectrum ( ), while in panels (b-d) the red line represents the target JONSWAP spectrum amended to water depth ℎ 2 by applying a shoaling coefficient from the linear wave theory.

Figure 6 :

 6 Figure 6: Case 1: spatial evolution of spectral parameters: (a) normalized peak frequency ̂ ∕ ; (b) spectral width parameter 1∕ ; and (c) normalized significant wave heights ∕ , for different spectral ranges. In all panels, the gray areas represent 95% CI. The vertical dash lines outline the extent of the bottom slope.

Figure 7 :

 7 Figure 7: Case 1: spatial evolution of statistical parameters: (a) skewness 3 , (b) asymmetry parameter and (c) kurtosis

Figure 8 :

 8 Figure 8: Evolution of the normalized wave spectra in cases 2-7, shown in panels (a-f) correspondingly. log 10 ( ( ) ∕ 0 ) is plotted in the colormap with a logarithmic scale. The close-to-horizontal dash line represents the evolution of ̂ and the vertical dash lines outline the extent of the bottom slope.

Figure 9 :

 9 Figure 9: Cases 2-7: Spatial evolution of the spectral parameters: (a) normalized peak frequency ̂ ∕ ; (b) spectral width parameter 1∕ ; (c) normalized significant wave height∕ , computed for different spectral ranges, defined in > 0.5 is shown in solid lines, < 0.5 in dash lines. The vertical dash lines outline the extent of the bottom slope, and the gray areas represent 95% CI.

Figure 10 :

 10 Figure 10: Cases 2-7: Equilibrium values of spectral parameters as functions of Ursell number ,1 of the incident sea-state, averaged in the region > 35 ,2 .

Figure 11 :

 11 Figure 11: Cases 2-7: Spatial evolution of the statistical parameters: (a) skewness 3 ; (b) asymmetry parameter ; (c) kurtosis 4 . The vertical dash lines outline the extent of the bottom slope, and the gray areas represent 95% CI.

Figure 12 :

 12 Figure 12: Cases 2-7: Maximum, minimum and equilibrium values of the statistical parameters as functions of Ursell number ,1 of the incident sea-state. The lines with stars represent the local extremes of the parameters achieved shortly after the slope, the lines with hollow circles represent the steady levels averaged in the region > 35 ,2 .

Figure 13 :Figure 14 :

 1314 Figure 13: Cases 2, 4 and 7: Evolution of cross-spectral parameters at characteristic positions in the short equilibration scale shown in upper panels. The corresponding coherence spectra are shown in the lower panels, where the horizontal dash line represents coh 2 1 , 2 = 0.25.

Table 1

 1 Summary of dimensional and non-dimensional wave parameters of the seven considered cases.

	Case	[m]	[m]	Deeper/shallower area
			1 ∕ 2	1 ∕ 2	,1 ∕ ,2
	1	0.030 103.6 1.85/0.64 0.037/0.062 0.0058/0.232
	2	0.010 43.6 1.85/0.64 0.012/0.020 0.0019/0.077
	3	0.015 43.6 1.85/0.64 0.019/0.031 0.0029/0.116
	4	0.020 43.

6 1.85/0.64 0.025/0.041 0.0039/0.155 5 0.025 43.6 1.85/0.64 0.031/0.052 0.0049/0.194 6 0.030 43.6 1.85/0.64 0.037/0.062 0.0058/0.232 7 0.035 43.6 1.85/0.64 0.043/0.073 0.0068/0.271

. In all panels the gray areas represent 95% CI, the vertical dash lines outline the extent of the bottom slope.

Acknowledgments

This work was supported by the China Postdoctoral Science Foundation (Grant No. 2021M690523), the National Natural Science Foundation of China (Grant No. 52101301), LiaoNing Revitalization Talents Program (Grant No. XLYC1807010) and the Fundamental Research Funds for the Central Universities (Grant No. DUT2019TB02).

Depth variations pose risks of freak waves only in the short scale. On the contrary, in the long scale, the depth variations could "prevent" structures or ships from freak waves to some extent. As freak wave occurrence probability decreases fast in the long scale, furthermore, high waves could be either filtered out by breaking effects or dissipated due to friction. A more dedicated investigation on the statistical distributions of free surface elevation, crest height and crestto-trough height during the equilibration process is presently ongoing and will be reported in an upcoming paper.

CRediT authorship contribution statement

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Methods of wave generation and absorption in the numerical model A.1 Wave generation

In the generation zone, the target incident wave field, denoted as ( ( , ), ̃ ( , )), can be calculated with either a linear or a nonlinear wave theory. Following the analysis of ZB21 and given the irregular wave trains of the cases listed in Table 1 are of relatively low level of nonlinearity in the deeper region, the incident wave field is generated here using the linear superposition technique. Considering a large number of ( ) harmonic components, the incident time series for wave-making is obtained with

where = 2 is the ℎ frequency component; is the corresponding wavenumber obtained via the dispersion relationship 2 = tanh ( ℎ 1 ); 0 and ℎ 1 correspond to the position and water depth at the wavemaker; is a random phase uniformly distributed in [0, 2 ); wave amplitude that corresponds to is chosen on the basis of the target energy density spectrum ( ), following eq. ( 1), defined in the range [ , ]. The target spectrum is divided into segments (here = 2 15 bins are set in [ , ]), each with the same portion of the total energy . Such a method renders unevenly distributed frequency components and the same wave amplitude for all components:

In the present study, we set = 0.45 and = 5 . During the numerical calculation, the actual values of the physical variables in the relaxation zone are updated at the end of each time step as

where stands for either or ̃ , ( ) is a ramp function progressively releasing the target solutions in space. For wave generation (left relaxation zone in Fig. 2), corresponds to the incident wave field. In addition, a Dirichlet condition of coefficients is imposed on the left lateral boundary of the computational domain, computed from the target wave solution. In the present application, the coefficient in the generation zone reads

A.2 Wave absorption

For the wave damping zone (right relaxation zone in Fig. 2), we use the same relaxation method as eq. (A.3) but now with = 0 to progressively attenuate the waves in space and limit wave reflection.

The coefficient is chosen as

Note that , varies from 0 at , = to 0.1 at , = + , in the damping zone. According to a number of calibration tests (not shown here), such a choice of , function in combination with a large , ends up with a good damping effect on long waves. In contrast to ZB21, where a different form of , ( ) was used and , was set to approximately 3 ,2 , we chose here a much longer zone with , = 21 m ≈ 20 ,2 in order to properly attenuate not only the wave components in the spectral peak range, but