Statistical dynamical model to predict extreme events and anomalous features in shallow water waves with abrupt depth change.

Understanding and predicting extreme events and their anomalous statistics in complex nonlinear systems are a grand challenge in climate, material, and neuroscience as well as for engineering design. Recent laboratory experiments in weakly turbulent shallow water reveal a remarkable transition from Gaussian to anomalous behavior as surface waves cross an abrupt depth change (ADC). Downstream of the ADC, probability density functions of surface displacement exhibit strong positive skewness accompanied by an elevated level of extreme events. Here, we develop a statistical dynamical model to explain and quantitatively predict the above anomalous statistical behavior as experimental control parameters are varied. The first step is to use incoming and outgoing truncated Korteweg-de Vries (TKdV) equations matched in time at the ADC. The TKdV equation is a Hamiltonian system, which induces incoming and outgoing statistical Gibbs invariant measures. The statistical matching of the known nearly Gaussian incoming Gibbs state at the ADC completely determines the predicted anomalous outgoing Gibbs state, which can be calculated by a simple sampling algorithm verified by direct numerical simulations, and successfully captures key features of the experiment. There is even an analytic formula for the anomalous outgoing skewness. The strategy here should be useful for predicting extreme anomalous statistical behavior in other dispersive media.

Understanding and predicting extreme events and their anomalous statistics in complex nonlinear systems are a grand challenge in climate, material, and neuroscience as well as for engineering design. This is a very active contemporary topic in applied mathematics with qualitative and quantitative models (1–7) and novel numerical algorithms, which overcome the curse of dimensionality for extreme event prediction in large complex systems (2, 8–11). The occurrence of Rogue waves as extreme events within different physical settings of deep water (12–16) and shallow water (17–19) is an important practical topic.

Recent laboratory experiments in weakly turbulent shallow water reveal a remarkable transition from Gaussian to anomalous behavior as surface waves cross an abrupt depth change (ADC). A normally distributed incoming wave train, upstream of the ADC, transitions to a highly non-Gaussian outgoing wave train, downstream of the ADC, that exhibits large positive skewness of the surface height and more frequent extreme events (20). Here, we develop a statistical dynamical model to explain and quantitatively predict this anomalous behavior as experimental control parameters are varied. The first step is to use incoming and outgoing truncated Korteweg–de Vries (TKdV) equations matched in time at the ADC. The TKdV equation is a Hamiltonian system, which induces incoming and outgoing Gibbs invariant measures. The statistical matching of the known nearly Gaussian incoming Gibbs state at the ADC completely determines the predicted anomalous outgoing Gibbs state, which can be calculated by a simple Markov chain Monte Carlo algorithm verified by direct numerical simulations, and successfully captures key features of the experiment. There is even an analytic formula for the anomalous outgoing skewness. The strategy here should be useful for predicting extreme anomalous statistical behavior in other dispersive media in different settings (21, 22).

1. Experiments Showing Anomalous Wave Statistics Induced by an ADC

Controlled laboratory experiments were carried out in ref. 20 to examine the statistical behavior of surface waves crossing an ADC. In these experiments, nearly unidirectional waves are generated by a paddle wheel and propagate through a long, narrow wave tank. Midway through, the waves encounter a step in the bottom topography and thus, abruptly transition to a shallower depth. The paddle wheel is forced with a pseudorandom signal intended to mimic a Gaussian random sea upstream of the ADC. In particular, the paddle angle is specified aswhere the weights are Gaussian in spectral space with peak frequency and the phases are uniformly distributed random variables. The peak frequency gives rise to a characteristic wavelength , which can be estimated from the dispersion relation. The energy injected into the system is determined by the angle amplitude , which is the main control parameter varied in ref. 20. Optical measurements of the free surface reveal a number of surprising statistical features.

Distinct statistics are found between the incoming and outgoing wave disturbances: incoming waves display near-Gaussian statistics, while outgoing waves skew strongly toward positive displacement.

The degree of non-Gaussianity present in the outgoing waves depends on the injected energy : larger energies generate stronger skewness in the surface displacement probability density functions (PDFs) and more extreme events.

Compared with the incoming wave train, the power spectrum of the outgoing wave field decays more slowly, which indicates that the anomalous behavior is associated with an elevated level of high frequencies.

2. Surface Wave Turbulence Modeled by TKdV Equation with Depth Dependence

The Korteweg–de Vries (KdV) equation is a one-dimensional, deterministic model capable of describing surface wave turbulence. More specifically, KdV is leading-order approximation for surface waves governed by a balance of nonlinear and dispersive effects valid in an appropriate far-field limit (23). Moreover, KdV has been adapted to describe waves propagating over variable depth (23). Here, we consider the variable-depth KdV equation truncated at wavenumber (with grid points) to generate weakly turbulent dynamics. The surface displacement is described by the state variable , with superscript − for the incoming waves and superscript for the outgoing waves. The Galerkin truncated variable is normalized with zero mean and unit energy , which are conserved quantities. Here, denotes the subspace projection. The evolution of is governed by the truncated KdV equation with depth change :Eq. is nondimensionalized on the periodic domain . The depth is assumed to be unit before the ADC and after the ADC. The conserved Hamiltonian can be decomposed aswhere we refer to as the cubic term and as the quadratic term. We introduce parameters based on the following assumptions.Some intuition for how Eq. produces different dynamics on either side of the ADC can be gained by considering the relative contributions of and in the Hamiltonian . The depth change, , increases the weight of and decreases that of , thereby promoting the effects of nonlinearity over dispersion and creating conditions favorable for extreme events. Since is the slope of the wave height, measures the wave slope energy.

The wavenumber truncation is fixed at a moderate value for generating weakly turbulent dynamics.

The state variable is normalized with zero mean and unit energy, , which are conserved during evolution. Meanwhile, characterizes the total energy injected into the system based on the driving amplitude .

The length scale of the system is chosen so that the resolved scale is comparable with the characteristic wave length from the experiments.

A deterministic matching condition is applied to the surface displacement to link the incoming and outgoing wave trains. Assuming that the ADC is met at , the matching condition is given byEq. is not designed to capture the short scale changes in rapid time. Rather, since we are interested in modeling statistics before and after the ADC, we will examine the long-time dynamics of large-scale structures.

Interpreting Experimental Parameters in the Dynamical Model.

The model parameters in [] can be directly linked to the basic scales from the physical problem. The important parameters that characterize the experiments (20) include , the wave amplitude to (upstream) water depth ratio; , the water depth to wavelength ratio; and , the depth ratio with upstream value and downstream value . By comparing the characteristic physical scales, the normalized TKdV Eq. can be linked to these experimental parameters viawhere defines the computational domain size as multiple of the characteristic wavelength , and , with a scaling factor for the state variable , normalizes the total energy of the system to 1.

Consider a spatial discretization with grid points so that the smallest resolved scale is comparable with the characteristic wavelengthIn the practical numerical simulations, we select and let vary in the range . Using reference values from the experiments (20), we find , and changes from 1 to 0.24 at the ADC. These values yield the following ranges for the model parameters: and . These are the values that we will test in the direct numerical simulations. Details about the derivation from scale analysis are in .

3. Equilibrium Statistical Mechanics for Generating the Stationary Invariant Measure

Since the TKdV equation satisfies the Liouville property, the equilibrium invariant measure can be described by a statistical formalism (24–26) based on a Gibbs measure with the conserved energy and Hamiltonian . The equilibrium invariant measure is dictated by the conservation laws in the TKdV equation. In the case of fixed energy , this is the mixed Gibbs measure with microcanonical energy and canonical Hamiltonian ensembles (24)where represents the “inverse temperature.” The distinct statistics in the upstream and downstream waves can be controlled with . As shown below, we find that negative temperature, , is the appropriate regime to predict the experiments. In the incoming wave field, is chosen so that has nearly Gaussian statistics. Using the above invariant measures [], the expectation of any functional can be computed asThe value of in the invariant measure is specified from (24, 26).

In addition to producing equilibrium PDFs of , the invariant measure can be used to predict the equilibrium energy spectrum without the direct simulation of TKdV. Direct simulation, however, is required to recover transient statistics of and time autocorrelations.

Statistical Matching Condition of the Invariant Measures Before and After the ADC.

The Gibbs measures are defined based on the different inverse temperatures on the two sides of the solutionswhere represents the deterministic matching condition between the incoming and outgoing waves. The two distributions, , should be matched at , givingIn matching the flow statistics before and after the ADC, we first use the conservation of the deterministic Hamiltonian after the depth change. Then, assuming ergodicity (24, 25), the statistical expectation for the Hamiltonian is conserved in time after the depth change at and should remain at this value as the system approaches equilibrium as . The final statistical matching condition to get the outgoing flow statistics with parameter can be found bywith the outgoing flow Hamiltonian and the Gibbs measures before and after the ADC.

4. Nearly Gaussian Incoming Statistical State

For the parameters explored in ref. 20, the incoming wave field is always characterized by a near-Gaussian distribution of the surface displacement. It is found that a physically consistent Gibbs measure should take negative values in the inverse temperature parameter , where a proper distribution function and a decaying energy spectrum are generated (ref. 26 and have explicit simulation results). The upstream Gibbs measure with displays a wide parameter regime in with near-Gaussian statistics. In Fig. 1, Left, the inflow skewness varies only slightly with changing values of and . The incoming flow PDF then can be determined by picking the proper parameter value in the near-Gaussian regime with small skewness. In contrast, the downstream Gibbs measure with shown in Fig. 1, Right generates much larger skewness as the absolute value of and the total energy level increases. The solid lines in Fig. 1 offer a further confirmation of the transition from near-Gaussian statistics with small to a strongly skewed distribution after the depth change.

Fig. 1.

(A) Skewness from the Gibbs measures in incoming and outgoing flow states with different values of total energy and inverse temperature (notice the different scales in the incoming and outgoing flows). (B) Outgoing flow parameter as a function of the incoming flow computed from the statistical matching condition with three energy levels . (C) Skewness in the outgoing flow with the matched value of as a function of the inflow parameter (the theoretical predictions using [] are compared).

In the next step, the value of the downstream is determined based on the matching condition []. The expectation about the incoming flow Gibbs measure can be calculated according to the predetermined parameter values of as well as from the previous step. For the direct numerical experiments shown later in Fig. 2, we pick proper choices of test parameter values as , and . More test cases with different system energy can be found in , where similar transitions from near-Gaussian symmetric PDFs to skewed PDFs in the flow state can always be observed.

Fig. 2.

Changes in the statistics of the flow state going through the ADC. The initial ensemble is set with the incoming flow Gibbs measure with different inverse temperature . (Row 1) Time evolution of the skewness and kurtosis. The ADC is taking place at . (Row 2) Inflow and outflow PDFs of . (Row 3) The downstream PDFs fitted with gamma distributions with consistent variance and skewness (in log coordinate in ). (Row 4) Energy spectra in the incoming and outgoing flows.

Direct Numerical Model Simulations.

Other than the prediction of equilibrium statistical measures from the equilibrium statistical approach, another way to predict the downstream model statistics is through running the dynamical model [] directly. The TKdV equation is found to be ergodic with proper mixing property as measured by the decay of autocorrelations as long as the system starts from a negative inverse temperature state as described before. For direct numerical simulations of the TKdV equations, a proper symplectic integrator is required to guarantee that the Hamiltonian and energy are conserved in time. It is crucial to use the symplectic scheme to guarantee the exact conservation of the energy and Hamiltonian, since they are playing the central role in generating the invariant measure and the statistical matching. The symplectic scheme used here for the time integration of the equation is the fourth-order midpoint method (27). Details about the mixing properties from different initial states and the numerical algorithm are described in .

5. Predicting Extreme Anomalous Behavior After the ADC by Statistical Matching

With the inflow statistics well described and the numerical scheme set up, we are able to predict the downstream anomalous statistics starting from the near-Gaussian incoming flow going through the ADC from to . We first consider the statistical prediction in the downstream equilibrium measure directly from the matching condition. The downstream parameter value is determined by solving the nonlinear Eq. as a function of , . In the numerical approach, we adopt a modified secant method, avoiding the stiffness in the parameter regime ( shows the algorithm). The fitted solution is plotted in Fig. 1 as a function of the proposed inflow . A nonlinear relation is discovered from the matching condition. The downstream inverse temperature will finally saturate at some level. The corresponding downstream skewness of the wave displacement predicted from the statistical matching of Gibbs measures is plotted in Fig. 1. In general, a large positive skewness for outgoing flow is predicted from the theory, while the incoming flow skewness is kept in a small value in a wide range of . Note that, with (that is, using the microcanonical ensemble only with energy conservation), the outflow statistics are also near Gaussian with weak skewness. The skewness in the outflow statistics grows as the inflow parameter value increases in amplitude.

For a second approach, we can use direct numerical simulations starting from the initial state sampled from the incoming flow Gibbs measure and check the transient changes in the model statistics. Fig. 2 illustrates the change of statistics as the flow goes through the ADC. Fig. 2, row 1 plots the changes in the skewness and kurtosis for the state variable after the depth change at . The PDFs in the incoming and outgoing flow states are compared with three different initial inverse temperatures . After the depth changes to abruptly at , both the skewness and kurtosis jump to a much larger value in a short time, implying the rapid transition to a highly skewed non-Gaussian statistical regime after the depth change. Furthermore, from Fig. 2, different initial skewness (but all relatively small) is set due to the various values of . With small , the change in the skewness is not very obvious (Fig. 2, row 2 shows the incoming and outgoing PDFs of ). In comparison, if the incoming flow starts from the initial parameter and , a much larger increase in the skewness is induced from the ADC. Furthermore, in the detailed plots in Fig. 2, row 3 for the downstream PDFs under logarithmic scale, fat tails toward the positive direction can be observed, which represent the extreme events in the downstream flow (Fig. 3 shows the time series of ).

Fig. 3.

Realization of the downstream (A) and upstream (B) flow solutions . Note the larger vertical scale in the downstream time series plot.

As a result, the downstream statistics in final equilibrium predicted from the direct numerical simulations here agree with the equilibrium statistical mechanics prediction illustrated in Fig. 1. The predictions from these two different approaches confirm each other.

6. Analytic Formula for the Upstream Skewness After the ADC

A statistical link between the upstream and downstream energy spectra can be found for an analytical prediction of the skewness in the flow state after the ADC. The skewness of the state variable at one spatial grid point is defined as the ratio between the third and second momentsNow, we introduce mild assumptions on the distribution functions.Then, the skewness of the downstream state variable after the ADC is given by the difference between the inflow and outflow wave slope energy of :The detailed derivation is shown in . In particular, the downstream skewness with near-Gaussian inflow statistics is positive if and only if the difference of the incoming and outgoing wave slope energy is positive. This means that there is more small-scale wave slope energy in the outgoing state. As evidence, in Fig. 2, row 4, in all of the weak and strong skewness cases, the outflow energy spectrum always has a slower decay rate than the inflow energy spectrum, which possesses stronger energy in larger scales and weaker energy in the smaller scales.

The upstream equilibrium measure has a relatively small skewness so that

The downstream equilibrium measure is homogeneous at each physical grid point so that the second and third moments are invariant at each grid point

In Fig. 1, we compare the accuracy of the theoretical estimation [] with numerical tests. In the regime with small incoming inverse temperature , the theoretical formula offers a quite accurate approximation of the third-order skewness using only information from the second-order moments of the wave-slope spectrum.

7. Key Features from Experiments Captured by the Statistical Dynamical Model

In this final section, we emphasize the crucial features generated by the statistical dynamical model [] by making comparison with the experimental observations in ref. 20. As from the scale analysis displayed in Section 2, the theory is set in the same parameter regime as the experimental setup.

The transition from near-Gaussian to skewed non-Gaussian distribution as well as the jump in both skewness and kurtosis observed in the experiment observations (figure 1 in ref. 20) can be characterized by the statistical model simulation results (Fig. 2, rows 1 and 2). Notice that the differences in the decay of third and fourth moments in the far end of the downstream regime from the experimental data are due to the dissipation effect in the flow from the wave absorbers that is not modeled in the statistical model here. The model simulation time series plotted in Fig. 3 can be compared with the observed time sequences from experiments (figure 1 in ref. 20). The downstream simulation generates waves with strong and frequent intermittency toward the positive displacement, while the upstream waves show symmetric displacements in two directions with, at most, small peaks in slow time. Even in the time series at a single location , the long-time variation displays similar structures.

The downstream PDFs in experimental data are estimated with a gamma distribution in figure 2 in ref. 20. Here, in the same way, we can fit the highly skewed outgoing flow PDFs from the numerical results with the gamma distribution

The parameters in the gamma distribution are fitted according to the measured statistics in skewness and variance: that is, . Additionally, the excess kurtosis of the gamma distribution can be recovered as . As shown in Fig. 2, row 3, excellent agreement in the PDFs with the gamma distributions is reached in consistency with the experimental data observations. The accuracy with this approximation increases as the initial inverse temperature increases in value to generate more skewed distribution functions.

Experimental measurements of the power spectra (figure 4 in ref. 20) reveal the downstream measurements to contain more energy at small scales (i.e., a relatively slower decay rate of the spectrum). This result is also observed in the direct numerical simulations here (detailed results are shown in ), as the outgoing state contains more energetic high frequencies.

8. Concluding Discussion

We have developed a statistical dynamical model to explain and predict extreme events and anomalous features of shallow water waves crossing an ADC. The theory is based on the dynamical modeling strategy consisting of the TKdV equation matched at the ADC with conservation of energy and Hamiltonian. Predictions can be made of the extreme events and anomalous features by matching incoming and outgoing statistical Gibbs measures before and after the abrupt depth transition. The statistical matching of the known nearly Gaussian incoming Gibbs state completely determines the predicted anomalous outgoing Gibbs state, which can be calculated by a simple sampling algorithm verified by direct numerical simulations, and successfully captures key features of the experiment. An analytic formula for the anomalous outgoing skewness is also derived. The strategy here should be useful for predicting extreme statistical events in other dispersive media in different settings.