Numerical modeling of rogue waves in coastal waters

Introduction Conclusions References


Introduction
This paper is dedicated to the analysis of rogue waves registered by buoys near the coast of Taiwan and to reproduction of these events with the help of hydrodynamic equations.The long-term measurements of the surface displacement by several buoys moored at different locations around Taiwan resulted in a collection of time series characterizing wave conditions in Taiwanese coastal waters (Lee et al., 2011).This data set was analyzed in Doong et al. (2012); in particular, rogue waves (the heights of which at least twice exceed the background significant wave height) have been identified and analyzed.The rogue (or freak) waves are a threat that has been recognized rather recently, and nowadays attracts much attention (see Kharif and Pelinovsky, 2003;Dysthe et al., 2008;Kharif et al., 2009;Slunyaev et al., 2011, for reviews).
Measurements of the surface elevation at a certain point provide only limited information about waves.The resulting time series, even if they contain rogue waves, do not reflect the dynamics of these dangerous waves.Much more information can be derived by merging the measured data with quite reasonable assumptions of wave dynamics.A common assumption is that waves propagate almost unidirectionally in a wind sea and particularly in swells.This implies that the angular spectrum is narrow.It is sufficient to fulfill this requirement within some reasonably short period of time, for example, for the duration of surface elevation record (10-20 min), or even shorter.Also, it is often true that the wave system contains no components (for example, reflected waves) propagating against the highest waves.Under these assumptions the surface wave dynamics can be considered as one-dimensional, for which the boundary value problem becomes in principle treatable.Numerical reconstruction of the evolution of such wave systems on the basis of instrumental records was performed in Trulsen (2001), Divinsky et al. (2004), Slunyaev et al. (2005Slunyaev et al. ( , 2014)), and Slunyaev (2006) by means of approximate evolution equations for modulated waves.As a result, not only can much more detailed wave information in the vicinity of the measurement point be obtained, but also the rogue wave evolution can be recovered.Importantly, application of the weakly nonlinear theory for wave modulations is able to produce reliable results even for very steep modulated waves (Slunyaev et al., 2014).A credible reconstruction of the wave dynamics was possible for about 10 min.
An important feature of the measured time series, which are in the focus of the present paper, is variable bathymetry: the sea is considered as a finite-depth basin, and the depth variation is essential.The complicated picture of changes in nonlinear waves when they travel from deep to shallow water was revealed by Sergeeva et al. (2011), Zeng andTrulsen (2012), andTrulsen et al. (2012).Though nonlinear wave interactions in constant-depth water are believed to trigger more rogue waves than predicted by the linear quasi-Gaussian statistics, these publications proved that variabledepth conditions were able to further increase the probability of rogue waves.
In this paper we consider four time series that are measured in the coastal waters of Taiwan and that contain rogue waves.In order to understand the nonlinear dynamics of extreme waves approaching from the open sea to the coast and to assess the rogue wave danger at the coast, we reconstruct the wave evolution of the in situ measured rogue waves employing the nonlinear theory for modulated waves.In Sect. 2 the conditions of measurements and the field data are briefly described.The employed evolution equations and the procedure of using the in situ time series to initialize the simulation are introduced in Sect.3. Section 4 provides the results of the simulations.The outcome is summarized in Discussion.

Field data
Long-term instrumental measurements of waves near the coasts of Taiwan have been conducted since 1997 by Taiwanese Central Weather Bureau (CWB) and Water Resources Agency (WRA) (Doong et al., 2007).Currently there are 15 buoys in operation in Taiwanese waters.The records used in this study are retrieved from four buoys moored at several locations, which differ in the distance from the coast, local depth, currents, typical wave spectra, etc.The sampling rate of the buoys is 2 Hz.The 512 s-long sections of time series of surface elevations are produced from the raw data of buoy accelerations with use of wavelet transform (see details in Lee et al., 2011).An extensive analysis of the available data set (Chang, 2012) exhibits a rich variety of rogue wave appearances.There are more than 100 rogue waves measured at each of those buoys.The occurrence probability of rogue  waves is in the range of 1.7-2.1 × 10 −4 ; that is, one rogue wave occurs in every 5000 waves.During the presence of a typhoon (tropical cyclone) close to the station, the frequency of occurrence of rogue waves increased by 2-7 times compared to its average level.It was also found that the percentage of rogue waves with symmetric shapes is large.
Four samples of the time series from the data set were picked out for the present study (see Figs. 1 and 2).We used the conventional criterion for individual rogue waves, where H is the height of a single wave and H s is the significant wave height.Note that we consider short segments of surface elevations (about 8.5 min instead of conventional 20 min), which may result in some overestimation of H s , and hence may lead to ignoring some rogue events.The quantity   AI will be referred to as abnormality index, which characterizes the deviation of the wave height from the significant value.The rogue waves represent three kinds of wave geometry: so-called "holes in the sea" with characteristic deep troughs (Fig. 2a, c, Hsinchu and Eluanbi stations), signvariable waves (Fig. 2b, Hualien station) and high-crested waves (Fig. 2d, Longdong station).The distances to shore, −x 0 , for the measuring stations are given in Table 1 together with the local depth h and other related parameters.Periods of the rogue waves are T fr = 6-10 s, where the subscript fr will denote the rogue (or freak) wave characteristics.The peak periods, T p , are obtained from the frequency spectra, whereas every record contains from 53 to 78 periods T p .Three buoys out of four are located at intermediate depths with dimensionless depths, k p h, within the range 1.3-2.3.The fourth buoy (Longdong) is moored at a deeper location with k p h = 4.Here k p is determined from the peak angular frequency, ω p , through the dispersion relation for finite-depth gravity waves, The significant wave heights, rogue wave heights and crest heights, H cr , are given in Table 1.The characteristic wave steepness, estimated from the standard deviation η rms of the surface elevation, as ε = k p η rms , varies from 0.02 to 0.04, which corresponds to rather weak nonlinearity.The steepness of the rogue wave may be estimated with the help of the local wave length λ fr (which is related to the local zerocrossing wave period according to the dispersion relation Eq. 2) and the crest height.The steepness of identified rogue waves k fr H cr is from 0.08 to 0.15.Therefore, the identified rogue waves were essentially weakly nonlinear.We used realistic bathymetry for numerical simulations of the wave propagation both onshore and offshore of the buoy locations.For each buoy the bathymetry was provided at a few locations (see blue circles in Fig. 3 and axis from the left side of the panels for the depth values) and was interpolated for our purposes using polynomial functions.

Simulation of weakly nonlinear waves over variable bottom
The time series of surface elevation, measured by individual gauges, contain very limited information.They basically present the local, instantaneous reaction of the water surface to the passing waves, but the wave dynamics remains concealed.To reconstruct the wave evolution in the vicinity of the measurement point, we employ weakly nonlinear equations that are able to describe the transformation of waves propagating in one direction (shoreward) over variable bathymetry.A suitable model is the nonlinear Schrödinger equation (NLS) for the case of wave evolution in space over uneven depth.The impact of bathymetry is taken into account through depth-dependent coefficients and a shoaling term (Djordjevic and Redekopp, 1978;Zeng and Trulsen, 2012): .
Here B(x, t) is the complex amplitude of the envelope, x is the horizontal coordinate directed onshore (x = 0 corresponds to the shoreline and the buoy is situated at x = x 0 ), t denotes the time, g is the gravity acceleration, and c p and c g are the local phase and group velocities, respectively.When the water depth h is constant, Eq. ( 3) reduces to the classic NLS equation for waves in a constant-depth basin.In this paper we assume that waves propagate strictly onshore with wave crests aligned with depth isolines.We follow the common concept that while waves propagate over variable depth, the carrier frequency ω 0 remains constant, and the wave number k varies in accordance with the dispersion relation Eq. ( 2); hence it is a function of x, k = k(h(x)).The dimensionless water depth kh is also a function of x (Fig. 3).As noted above, the local dimensionless depth corresponds either to the intermediate or to relatively deep water at x 0 , and corresponds to the truly deep water at the farthest modeled position offshore.The ocean at the Hsinchu station is shallower.
Note that when ω 0 is constant, the group velocity c g depends on the x coordinate through the product k(x)h(x); see expression for c g in Eq. (3).Therefore, the coefficient µ may be written in the form The shoaling term (the first additive at the right-hand side of Eq. 3) may be expressed as −i/(2c g ) dc g /dx; hence it reflects the conservation of wave energy flux c g |B| 2 as could be expected from the linear theory.
The numerical code for solving the NLS Eq. ( 3) is implemented similarly to Lo and Mei (1985), Slunyaev et al. (2005), Slunyaev (2006), and Shemer et al. (2010).We use the split-step-Fourier method, and the integration in space is performed for the periodic in time domain using the fourth-order Runge-Kutta method.The classical integration in space represents the forward wave evolution, whereas the backward evolution is obtained using a negative step of integration in space.
The measured time series are used as the starting condition.The carrier frequency is set equal to the mean frequency of the wave energy spectrum.The complex envelope B(x 0 ), at x = x 0 , is sought with the help of an iterative procedure.This procedure assumes that the in situ surface  for weakly nonlinear weakly modulated wave trains.This assumption means that the second-and third-order wave corrections and also the large-scale wave-induced flow are taken into account (see details of the approach in Trulsen (2001Trulsen ( , 2006)), Slunyaev (2005), and Slunyaev et al. (2014)).Similarly, the bound wave harmonics are computed for a given B(t) to produce the surface elevation function η(t).Strictly speaking, for relatively shallow water (kh 2), coefficients in the reconstruction formulas differ from those in the deepwater limit.For intermediate depth the coefficients are extremely cumbersome (Slunyaev, 2005); that is why a simplified approach was employed in this study.

Reconstruction of rogue waves by means of the variable-coefficient NLS equation
The essence of our simulations is that the wave evolution in the onshore direction (increasing x) and offshore from the buoy location is reconstructed from the measured time series within the framework of the variable-coefficient NLS as described above.When waves approach the shoreline, their amplitude grows, and at some location the wave steepness becomes unrealistically large.The wave breaking is not resolved by the weakly nonlinear model, and waves can grow infinitively as h → 0 in this model.For this reason simulations are only performed until the location of the shallowest available bathymetry data point (the leftmost circles in Fig. 3) -for example, until about 130 m from the coast for the Eluanbi buoy.The simulation recovers the wave evolu-tion backward for a few kilometers offshore of the buoys (Fig. 3).The simulated data were stored each 2 m of the wave run; the output data corresponding to wave fields of the size 512 s cover a spatial domain from 3.5 to 6 km.The rogue wave evolution is restored by virtue of the approach.These data may be used for a study of the rogue wave evolution and appearance, and of mutual properties of the elevation and kinematic fields (Sergeeva and Slunyaev, 2013).Since the local wave energy varies due to the shoaling effect, the significant wave height H s = H s (x) is calculated from the reconstructed time series at particular locations.The significant wave height is estimated in the quasi-Gaussian approximation from the standard deviation η rms of the surface elevation as H s = 4η rms (Massel, 1996;Kharif et al., 2009).The reconstructed rogue wave disappears at some distance from x 0 in the sense that its height falls below the level of 2H s (x).However, it may again exceed this threshold shortly after, and thus it is natural to consider the reappearing extremely large waves as one continuous rogue wave event.The lifetimes of such events, when the condition Eq. ( 1) may be violated for some short periods of time (less than a couple of wave periods), may reach up to 10 min in deep-sea simulations of unidirectional Euler equations (Sergeeva and Slunyaev 2013).
Such "intermittent" rogue events are observed in the present simulations.An example is provided in Fig. 4 for the record from the Hsinchu buoy.The rogue wave emerges a little bit before the buoy location (x 0 = −2440 m) and propagates for about 550 m for about 60 sec, from time to time exceeding the threshold 2H s .Quite naturally, the rogue wave shape in the simulated wave profiles changes in the course of evolution.Similar pictures are obtained for other recorded time series.Although time series in Fig. 4 demonstrate that the rogue waves often belong to intense wave groups, no cases were observed when subsequent waves in the time series exceeded the value of 2H s (i.e., no rogue groups).
The spatio-temporal (x, t) fields of the surface elevation (Fig. 5) reveal short-living, long-living and intermittent rogue waves (in the two latter cases waves are collated by red boxes).Note the periodic boundary condition for the time axis; thus when a wave reaches the right boundary of the shown domain, it continues from the left side.It is readily seen that most of the long-living rogue events occur on the background of intense trains that do not exhibit a clear dispersion of energy; this confirms observations in Sergeeva and Slunyaev (2013) for the deep-sea conditions.
Eight rogue events satisfying Eq. (1) were found in the simulations (Fig. 5).Four of them correspond to the rogue waves registered at x 0 , and the others occur at other locations.Importantly, all the "new" rogue events occur seaward of the buoys (i.e., at deeper water), whereas no "new" rogue waves were detected in the shallower areas.The rogue event lifetimes and corresponding distances passed by the rogue structures are given in Fig. 6 as functions of the abnormality index, AI= H /H s .There are long-living extreme events that last up to 100 s and travel 600 m, and also transient events.Generally, rogue waves with large AI live longer.
Each simulated rogue event contains up to 300 time series enclosing rogue waves.This representation makes it possible to estimate the most probable appearance of the rogue wave in statistical sense.Typical shapes of observed in situ rogue waves were discussed, for example, in Kharif et al. (2009) and Didenkulova (2011).Four kinds of rogue waves are distinguished: rear positive waves (RPWs), rear negative waves (RNWs), front positive waves (FPWs), front negative waves (FNWs).This classification was introduced in Sergeeva and Slunyaev ( 2013) and reflects the wave asymmetry.For example, the FPW has an extreme front slope and its crest height exceeds the depth of its trough.The RNW has the extreme slope at the rear side, and the trough is deeper than the height of the crest.
Figure 7 represents the frequency of occurrence of these four kinds of waves for each of the eight rogue events.Although the numbers are not always sufficiently large to provide a statistically sensible result, Fig. 7 suggests that the major part of rogue waves have high crests and shallow troughs.Waves with deeper troughs are observed as well (similar to shown in Fig. 1a), though they are less numerous.The proportion of waves with the extreme slope at the front or at the rear side of waves is about the same.

Discussion
The rogue wave evolution is reconstructed on the basis of in situ instrumental records of the surface elevation in the manner similar to Trulsen (2001), Slunyaev et al. (2005Slunyaev et al. ( , 2014)), and Slunyaev (2006).The asymptotic equations for surface water waves, which describe wave evolution in space, are solved numerically.Waves are assumed to be unidirectional, weakly modulated and weakly nonlinear.This approach is able to adequately replicate even fairly steep modulated waves (Slunyaev et al., 2014).The new aspect in our study is wave propagation over essentially variable bottom profile; hence we use a modification of the evolution equation that accounts for shoaling effects.
The recent numerical and partly laboratory studies demonstrated that the wave nonlinearity was able to increase significantly the probability of rogue waves both in shallow (Pelinovsky and Sergeeva, 2006) and deep water (Onorato et al., 2001(Onorato et al., , 2002;;Janssen, 2003;Shemer and Sergeeva, 2009, and   many subsequent publications).Hence, the nonlinear wave interactions in homogeneous conditions may lead to a higher risk of rogue wave occurrence.
When waves travel from deeper to shallower areas, the wave amplitude and energy is amplified due to the change of propagation velocity.This classical quasi-linear effect is captured by the NLS Eq. (3). Figure 8 shows the changes to the wave amplitude and energy when waves are approaching the coast, through the quantity 8η rms (thick black line) near the Hsinchu station.As the energy grows further when the depth decreases, at some location the numerical integration becomes inaccurate.
At first glance, it is reasonable to expect that this intensification of wave amplitude may result in a greater number of rogue waves when waves propagate shoreward.Numerical and laboratory simulations by Sergeeva et al. (2011), Zeng andTrulsen (2012), andTrulsen et al. (2012) of irregular nonlinear waves crossing a smooth step draw a more complicated picture, still to be comprehended.The variation in the water depth can lead to significant changes in statistical properties of the wave field, and relatively fast changes in the properties of the waveguide may lead to radical alteration of the wave dynamics.In particular, an increase in the statistical moments in the shallower region may enhance the likelihood of extreme wave heights.In some situations a significant increase in the rogue wave occurrence was observed indeed, though in other cases no significant changes were found.Zeng and Trulsen (2012) considered a "deeper" regime (kh > 1.2) of wave evolution compared to numerical simulations by Sergeeva et al. (2011) (kh < 0.4) and exposed the reduced kurtosis and skewness for the shallower region of transformation.Laboratory experiments of Trulsen et al. (2012) embrace the "transitional" zone when the waves travel from an intermediate depth to the shallow water with kh down to 0.54.An increase of large wave likelihood and kurtosis was observed for the smallest dimensionless depth kh = 0.54 in the experiment.The shallow water conditions of kh equal to 0.7 and 0.99 demonstrated the decrease of these statistical parameters.Let us reiterate an interesting peculiarity of gravity waves: their group velocity possesses a maximum at kh ≈ 1.2 and is somewhat smaller over deeper seas and much smaller in coastal waters.Thus, the linear theory predicts that waves propagating onshore from the deep water should somewhat diminish at kh ≈ 1.2.
When waves travel from deeper to shallower water, the problem is controlled by the dimensionless initial and final local depths, the Ursell parameter and the Benjamin-Feir index (BFI).The Benjamin-Feir index may be written in a form that explicitly reflects the finiteness of the water depth: Here ν and λ are coefficients specified in Eq. ( 3), and the rightmost factor in Eq. ( 5) tends to one in the limit of infinite depth.The BFI characterizes the importance of nonlinear effects due to four-wave resonances with respect to the wave dispersion.Uniform waves become modulationally unstable when BFI > 1.This index becomes less than zero when kh becomes smaller than 1.363 due to the change of sign of the coefficient ν (Johnson, 1997;Slunyaev et al., 2002).We emphasize that when the sea is deep (kh 1), envelope solitons are exact solutions of the NLS equation with constant coefficients.These solutions are wave groups that in the theory may propagate eternally, and for a rather long time in realistic conditions (see Slunyaev et al., 2013, and references therein).Long-living intense wave groups in irregular unidirectional waves over deep water are believed to be related to the NLS solitons (Slunyaev and Sergeeva, 2011).The long-living wave groups clearly observed in Fig. 6 resemble soliton-like structures that are transforming when propagating over variable depth.At the same time, the BFIs calculated for the typical wave conditions at the measurement locations are quite small (Table 1).When computing BFI, the spectral width is determined through the second spectral moment within the range 0.4-1.2rad s −1 (it corresponds approximately to the free wave spectral domain).The in situ wave spectra have different shapes, and generally are not narrow; hence generalized equations for wave modulations should be more accurate.
When the water depth decreases, the nonlinear coefficient in the NLS equation also decreases so that BFI decays as well.This process reduces the self-focusing effect.As a result, solitary groups of a given amplitude consist of a greater number of waves (see Slunyaev et al., 2002;Didenkulova et al., 2013).In even more shallower sea with kh < 1.363, solitary wave groups are no more solutions for the NLS equation.Wave conditions and water depth at the Hsinchu station correspond to kh < 1.363, and thus the Benjamin-Feir index has a negative value.As envelope solitons disperse when the depth becomes h < 1.363/k (Grimshaw and Annenkov, 2011), the role of solitary groups is expected to diminish as kh 2, when the coefficients of the NLS equation differ significantly from their values in deep water.The values of BFI observed in all the reported simulations remain less than 0.17.
As water becomes shallower, three-wave resonances start to play the major role.The importance of the shallow-water nonlinearity is controlled by the Ursell parameter: The Stokes wave theory requires Ur < 10 (see Holthuijsen, 2007).When Ur is too large, the NLS theory for modulated waves becomes inapplicable (when Ur > 26, according to Holthuijsen, 2007).This parameter is calculated for the measured time series taking in Eq. ( 6) H = 4η rms and ω = ω p .The values of Ur at the measurement locations are sufficiently small to employ the NLS theory (Table 1).In the shallow water limit kh 1 the conservation of energy flux results in Green's law in the form H (x) /H (x 0 ) = (h(x 0 )/h(x)) 1/4 , and the Ursell number grows as the sea becomes shallower.In the numerical simulations, the Ursell numbers remain very small (Ur < 2) in cases of Hualien, Eluanbi and Longdong buoys; in contrast, the waves from Hsinchu buoy correspond to a very large value Ur ≈ 38 at the shallowest available bathymetry location kh ≈ 0.66.All four numerical simulations reported by us exhibit the situation when coastal rogue waves occur more frequently at deeper locations (Fig. 5). Figure 8 presents the dependence of maximum wave height in the simulations of the variablecoefficient NLS equation (the thin black curve) versus dimensionless depth, kh.The bold line shows the dependence of 8η rms , which characterizes the evolution of the mean wave amplitude.The sea is rather shallow in the location of the Hsinchu buoy (k p h = 1.34).Despite the confident growth of the average wave height close to the coast, the wave heights do not exceed the limit 8η rms .
Due to the large number of controlling parameters, the conclusion that is drawn from the present simulations cannot be considered general: other combinations of wave intensity and bathymetry peculiarities might result in opposite effects; in particular, it is not excluded that in accordance with observations collected by Nikolkina and Didenkulova (2012) the rogue wave probability may increase as waves propagate shoreward.It is crucial that other physical effects may start to act as waves approach shallower regions -for example, interactions of shallow water solitary waves (Peterson et al., 2003).This effect is totally neglected by the theory that is employed in this paper.

Fig. 1 .
Fig. 1.Locations of the buoys used in this study.

Figure 3 .
Figure 3.The uneven bathymetry in the vicinity of the measurement locations.The circles and the connecting blue lines show the depth.The local dimensionless water depth kh is shown with broken red lines and corresponding right-side axis.The vertical lines denote the measurement positions x 0 .Locations: a) Hsinchu, b) Hualien, c) Elaunbi, d) Longdong.

Fig. 3 .
Fig. 3.The uneven bathymetry in the vicinity of the measurement locations.The circles and the connecting blue lines show the depth.The local dimensionless water depth kh is shown with broken red lines and corresponding right-side axis.The vertical lines denote the measurement positions x 0 .Locations: (a) Hsinchu, (b) Hualien, (c) Elaunbi, and (d) Longdong.

Figure 4 .
Figure 4.The simulated rogue wave evolution in the vicinity of the Hsinchu buoy.Numbers at the left side give values of x in meters.The thick red curves show rogue wave according to Eq. (1) for local x.

Fig. 4 .
Fig. 4. The simulated rogue wave evolution in the vicinity of the Hsinchu buoy.Numbers at the left side give values of x in meters.The thick red curves show rogue wave according to Eq. (1) for local x.

Fig. 5 .
Fig. 5.The spatio-temporal fields of surface elevations with rogue waves (blue circles), and long-living and intermittent rogue events (red boxes).Simulations of records from the Hsinchu (a), Hualien (b), Longdong (c) and Eluanbi (d) buoys.Horizontal lines indicate coordinates of the in situ measurements.

Figure 6 .
Figure 6.Life-times (a) and passed distances (b) of rogue events in the simulated data versus the abnormality index, AI.

Fig. 6 .Fig. 7 .
Fig. 6.Lifetimes (a) and passed distances (b) of rogue events in the simulated data versus the abnormality index, AI.

Fig. 8 .
Fig. 8. Maximum wave heights in meters versus dimensionless sea depth in the numerical simulations of the NLS; the thick lines give values of 8η rms .Simulations of record from the Hsinchu buoy are shown.

Table 1 .
Main parameters of the in situ records.

fr [m] H fr cr [m] H s [m] Ur
p [s] T fr [s]  fr [m] H