The ocean surface and, by extension, ocean wave energy are probabilistic and should be understood via probabilistic analysis. In the present work, which represents a seed that establishes a solid theoretical foundation on which the future work can be built, we demonstrate a probabilistic approach to the time evolution of ocean wave energy via a semi-analytic solution using the Wiener chaos expansion method. We present a comparison between field observations and corresponding Wiener chaos expansion calculations of the potential and kinetic energies of ocean surface waves. We also compare Wiener chaos expansion calculations of ocean surface kurtosis with field observations. Significant characteristics of the behavior of field-data are seen in the results produced by the Wiener chaos expansion method. This demonstrates the possibility of the use of the Wiener chaos expansion method in understanding the probabilistic behavior of the time-evolution of total ocean wave energy for capture by wave power devices.

## I. INTRODUCTION

Robust laboratory experiments^{1–7} and compelling field data evidence^{8–11} present indisputable evidence of non-Gaussianity in the ocean surface. The effect of this non-Gaussianity on ocean wave energy is our primary concern. It is widely accepted that ocean surface elevation statistics deviate from Gaussianity,^{10,12,13} and this is seen in numerous simulations of the ocean surface wave envelope^{14,15} as well as statistical approaches to the ocean surface,^{16–18} but the analytical connection between the direct ocean surface elevation momentum equations and the non-Gaussian behavior seen in wave height statistics, such as rogue wave occurrence, is the subject of ongoing research. Oceanic non-Gaussianity (where ocean surface elevation shows statistical distributions that vary significantly from the Gaussian distribution) has been associated with physical phenomena such as instability in wind-waves,^{12,19–23} nonlinear energy transfer between waves,^{12,19,21,24,25} and diffraction of near-unidirectional waves.^{15,25–29} We study non-Gaussianity resulting from the nonlinear time evolution of Gaussian basis functions in the direct moment equations for surface water waves. We then apply this method toward ocean surface energy calculations.

Spectral domain models, though computationally expensive, are more efficient than time domain models for simulating ocean waves, as the author previously demonstrates,^{30} where realistic ocean wave fields are simulated via higher order spectral methods using a Monte Carlo simulation approach. Designers who model in the spectral domain use probabilistic approaches such as Monte Carlo methods to model input waves as a Gaussian process with known statistical characteristics^{31,32} including that surface elevation is perfectly mean-centered with zero skewness (third statistical moment) and zero excess kurtosis (fourth statistical moment).^{33}

However, evidence of oceanic non-Gaussianity is usually found in the waves that result from these spectral approaches to ocean surface simulation, showing that a process that begins as Gaussian can produce non-Gaussian results based on the non-linearity of oceanic dynamics.

In the present study, a less computationally expensive stochastic approach to calculating the time evolution of the total potential and kinetic energies of the ocean surface is examined. Our primary (and novel) method of calculating stochastically driven ocean wave energy is via the use of Wiener Chaos Expansion (WCE). The WCE method we employ was established in its present form by Ghanem^{34} after being previously developed to solve probabilistic problems of applied mechanics by Cameron and Martin.^{35} In the WCE method, we represent ocean surface waves as expansions of Gaussian random variables. We directly solve for the WCE coefficients of ocean surface elevation via the Craig and Sulem^{36} formulation of the Zakharov equations, analytically arriving at calculations for ocean surface statistical moments and estimates of ocean wave energy.

The application of the present method is the capture of the non-Gaussian evolution of higher order statistical moments of the ocean surface via only one time-domain, deterministic system of governing equations. Herein, the kinetic energy and potential energy of the ocean surface are mathematically treated as higher order moments of the ocean surface. Thus, estimates of the time-domain variability of the total kinetic energy and total potential energy in an oceanic wave field are deterministically achieved. This can serve as an indicator of the total environmental energy that is available to wave power devices for capture.

Typically, other time-domain calculations require several simulations estimating average power capture and its standard deviation.^{31} We recommend the Wiener chaos expansion approach as an alternative source of environmental forcing, since it involves, as demonstrated below, only one simulation to also capture the temporal result evolution of the total usable mechanical energy in a wave field.

## II. WCE-BASED DEFINITION OF OCEAN SURFACE ELEVATION

If we denote the velocity potential by *ϕ*, where velocity **u** = ∇*ϕ*, then the velocity potential defined within the fluid domain (see Fig. 1) is governed by the continuity equation in the fluid domain,^{37}

where the free surface is denoted by *y* = *η*(*x*, *t*), the **y**-value of the mean of the fluid surface is set as *y* = 0, and the **y**-value of the bottom of the fluid domain is *y* = −*h*.

We enforce a no-flux boundary condition at the bottom, where *y* = −*h*, preventing vertical velocity at the bottom of the fluid domain,

For a gently sloping free surface, the vertical velocity of the fluid on the free surface equates to the vertical velocity of the surface itself, and the surface kinematic boundary condition is, thus, defined as the absence of a jump in velocity at the water surface,^{47}

We represent the surface dynamic boundary condition as the absence of a jump in pressure at the water surface, so that the pressure due to the local fluctuations in velocity potential and the hydrostatic pressure are balanced by the pressures due to convective effects in the fluid domain and static atmospheric pressure (zero-value on the left-hand side) as follows:

The intended stochastic analysis requires variables that are expressed in two dimensions. We, therefore, use the notation and formulation of Craig and Sulem,^{36} in which the water wave problem is described in terms of two surface-defined variables, namely, *η*, surface elevation, and *ξ* that represents the velocity potential defined at the free surface.

The Dirichlet-to-Neumann operator *G*(*η*) is an analytic function of *η* and can be defined as follows:

where *η* represents the free surface elevation and *ξ* represents the velocity potential, *ϕ*, at the free surface.^{36} With this definition of the Dirichlet-to-Neumann operator, *G*(*η*), Craig and Sulem describe the surface water wave equations in variables defined only at the water surface. Written according to the format of Craig and Sulem,^{36} we express the free-response water-surface wave equations as follows:

which is the kinematic boundary condition, and

which is the free-response dynamic boundary condition.

Given the nonlinear terms of Eq. (7), a truncated Taylor expansion is used to express the $121+\eta x2$ coefficient as $121\u2212\eta x2+\eta x22\u2212\cdots $, and this is placed in the free-response dynamic boundary condition as follows:

dropping the negligible higher order terms and keeping only the underlined terms in the coefficient term.

Furthermore, to represent surface waves driven by random wind, we add a random forcing term to the free-response dynamic boundary condition as done by Craig and Sulem^{36} in their ship-load example. The main governing equations are then re-written as follows:

which remains as the kinematic boundary condition, and

which becomes the randomly driven dynamic boundary condition, where *σ* represents the spatial distribution of random forcing and *W*(*t*) represents a Brownian motion forcing term.

The Dirichlet-to-Neumann operator, *G*(*η*), can be described in the form of a convergent Taylor expansion, whose first term is defined as

where the Fourier transform of the convolution of a variable with −*i∂*_{x} is mathematically equivalent to pre-multiplication by *k* in Fourier space.^{38,39}

The second and third terms are

and

respectively.

We use only the first term *G*_{0} in our estimations. This is justified in Figs. 2–4 where we simulate a second order Stokes wave with initial conditions

and

where

and

Note also that wavenumber *k* = 3 rad/m, frequency $\omega =gktanh(kh)$, fluid domain depth *h* = 1 m, acceleration due to gravity *g* = 9.81 m/s^{2}, and dimensionless constants *β* = 0.74 and *α* = 8.1 * 10^{−3}.

These results demonstrate that higher order terms in the Dirichlet-to-Neumann operator *G*(*η*) do not significantly modify the ocean surface and kurtosis estimations.

Our numerical tests confirm the adequacy of using only the first term *G*_{0} as an approximation of the function of *G*(*η*). We, therefore, use only the first term *G*_{0} in our present estimations.

Randomly driven surface variables *η* and *ξ* expressed in this two-dimensional form can, thus, be deduced via an infinite Fourier–Hermite Wiener chaos expansion system. *η* and *ξ* can, therefore, be expanded as follows:

via infinite expansions that are truncated for computational tractability. As shown in Eqs. (19) and (20), each term in the infinite system is labeled with a specified index, *α*. These indices are infinite in the infinite system but are limited in the truncated system. Truncation order is denoted by two numbers:

*M*, the order of nonlinearity of the system [equivalent to the maximum allowed sum of the individual sub-index terms*α*_{1},*α*_{2}, …,*α*_{K}among the allowed indices*α*= (*α*_{1},*α*_{2}, …,*α*_{K})], and*K*, the number of random kicks allowed in the truncation (equivalent to the number of sub-index terms in each index).

We can denote the entire set of applicable indices *α* = (*α*_{1}, *α*_{2}, …, *α*_{K}) in a specified truncation as $J$. For each applicable index in a truncated problem, the relevant multi-variable Hermite polynomial (or Wick polynomial) *ζ* = (*ζ*_{1}, *ζ*_{2}, …, *ζ*_{K}) is defined as

where we define the independent identically distributed functions *ζ*_{i} as random variables distributed according to the standard normal distribution. The $H\alpha i(\zeta i)$ term in Eq. (21) is a Hermite polynomial. The generating function for a Hermite polynomial is defined via Taylor series as

The Wiener chaos expansions of surface elevation, *η*(*x*, *t*), and surface velocity potential, *ξ*(*x*, *t*), have the following Fourier–Hermite expansions:

where the coefficients are the *η*_{α}(*x*, *t*) and *ξ*_{α}(*x*, *t*) terms.

Similarly to the work of Hou,^{39} we can solve for the time evolution of this system’s Wiener chaos expansion coefficients by applying them to the above stated nonlinear ocean surface elevation [Eqs. (9) and (10)] as follows:

where *p*, *β*, and *α* are allowable multi-indices in the doubly truncated Wiener chaos expansion of the ocean surface elevation equations and *C*(*α*, *β*, *p*), the total number of permissible combinations within the truncation, is deduced via

*α*, *β*, *p*, *γ*, and *q* are all *K*-long, *M* (maximum sum)-multi-indices in a chosen *M*-*K* truncation, such that

*α*≥*β*,$p\u2208J$,

$q\u2208J$,

*α*−*β*+*p*≥*γ*,$\alpha +2p\u2208J$, and

$\alpha \u2212\beta +p+2q\u2208J$.

The deterministic WCE representation of the random-forcing term, $\sigma (x)W\u0307(t)$ in Eq. (25b), is represented by $\sigma (x)\u2211i=0\u221eI\alpha i=\delta ijmi(t)$. For fixed *i*,

We define the orthogonal basis functions of Brownian motion effects in the water wave WCE *m*_{k}(*s*) as follows:

where 0 ≤ *s* ≤ *t* and 2 ≤ *κ* ≤ *K*, in which *K* depends on the order of truncation of the Fourier–Hermite system. Together,

represent the deterministic effects of random surface wind as the first-order terms in the truncated WCE system (see further details in Ref. 39).

## III. WCE-BASED DEFINITION OF OCEAN SURFACE KURTOSIS

Kurtosis is commonly referred to as the fatness of the tails in a probability distribution function. It is defined as a measure of the likelihood of extreme values in a distribution. The kurtosis *κ*_{4} of a theoretical distribution is expressed in the following relationship between the fourth statistical moment *μ*_{4} and the second statistical moment *μ*_{2}:

The kurtosis of a normal (Gaussian) distribution is equal to 3. “Excess kurtosis” is, therefore, often defined as *κ*_{40} = *κ*_{4} − 3 and is considered a standard indicator of non-Gaussianity. When a distribution has a value of *κ*_{40} < 0, it is considered platykurtic. Comparably, when it has a value of *κ*_{40} = 0, it is considered mesokurtic, and if *κ*_{40} > 0, it is considered leptokurtic.

Many authors observe ocean surface elevation kurtosis via the Zakharov equations.^{19,39–44} The Zakharov equation is a differential equation governing the behavior of a canonical complex variable, **b**(**k**), whose real and imaginary parts represent the Fourier transforms of surface elevation, $\eta \u0302k,t$, and of surface velocity potential, $\xi \u0302k,t$, respectively.

Solutions of Zakharov’s equation prove 2D, 3D, resonant, non-resonant, random, and non-random scenarios in which changes in kurtosis arise.

In the present work, using values calculated from the solution of Eqs. (25a) and (25b), ocean surface kurtosis is defined in terms of Wiener chaos expansion coefficients of the ocean surface *η* as follows:

## IV. WCE-BASED DEFINITION OF OCEAN SURFACE POTENTIAL ENERGY

According to linear theory, the energy density, E, in an ocean wave field of surface elevation *η* in depth h is the sum of its kinetic energy and potential energy^{37} (p. 260). In a near-linear wave field, ignoring energy losses, this energy density is best approximated as evenly split between the potential energy and the kinetic energy. For a homogeneous wave field in deep water, the mean potential energy per unit area^{45} (p. 521), *V*, can be expressed as

Following this, we use the WCE method^{34} according to the convention of Hou *et al.*^{39} to express the Wiener chaos expansions of ocean surface elevation and ocean surface velocity potential. We express the WCE coefficients of the potential energy, *V*_{α}, in terms of the calculated WCE coefficients of ocean surface elevation, *η*_{α},

The ocean surface elevation WCE potential energy coefficients *V*_{α}’s can then be summed as follows:

where *T*_{α} is the *α*-order Wick polynomial. This is used to calculate an analytical representation of ocean surface elevation potential energy, V, as a function of time.

Starting with a second-order Stokes wave as the initial condition, we apply the Wiener chaos expansion at truncation order K = 3, M = 2. Estimations for *η*_{α} are used to calculate the time evolution of potential energy, V. The results are shown in Fig. 5.

## V. WCE-BASED DEFINITION OF KINETIC ENERGY

Similarly, the mean kinetic energy (*KE*) of a homogeneous wave field in deep water can be defined as^{45} (p. 522)

where the surface-defined velocity potential of the fluid domain *ξ* = *ϕ*(*x*, *y* = *η*, *t*) and the surface elevation *η* = *η*(*x*, *t*).

This leads to the WCE expression for the coefficient of kinetic energy, *KE*, in terms of the calculated WCE coefficients *η*_{α} and *ξ*_{α},

We can then approximate that the total kinetic energy *KE* is as follows:

where *T*_{α} is a Wick polynomial as defined and described in Sec. V.

Starting with a second-order Stokes wave as the initial condition, we apply the Wiener chaos expansion on the Stokes wave at truncation order K = 3, M = 2. Estimations for *η*_{α} and *ξ*_{α} are used to calculate the time evolution of kinetic energy, *KE*. The results are shown in Fig. 6.

## VI. DATA COLLECTION FOR WCE METHOD INITIAL CONDITION

In the Las Cuevas Bay, off the North Coast of Trinidad, a 1200 kHz Teledyne Acoustic Doppler Current Profiler (ADCP) Workhorse Sentinel device is deployed in the seabed at 14.3 m water depth, as depicted in Fig. 7. The ADCP is deployed upward-facing and lodged to the ocean floor with a heavily weighted fiberglass anchor.

The ADCP measures surface elevations and surface velocities via measurements of the echo-Doppler effect along beams pinged from the device transmitters (indicated via red ovals in Fig. 7) to the surface and measured on their way back as they reflect off the ocean surface back to the device receptors (also indicated via red ovals in Fig. 7). For the relevant frequency range of waves, seven hours of 2 Hz data allow for a full understanding of the pertinent spectral frequencies.

Data clean-up is done according to the method of Boyd,^{46} taking the following steps:

Remove erroneous values that are found in small amounts in each dataset,

Replace erroneous values with an interpolated value being the average of the adjacent non-zero-values,

High-pass filter data to remove the tidal effect.

Starting with a 20-minute record of Las Cuevas Bay surface elevation and surface velocity data, we define the initial conditions of the lowest order terms in the Wiener chaos expansions.

The lowest order coefficient of surface elevation $\eta 000*$ is defined with the following dimensionless term:

where *η* is defined by surface elevation data.

For the lowest order velocity potential terms, we make substitutions of terms at the relevant orders into Eq. (32). The lowest order coefficient of the surface-defined velocity potential $\xi 000*$ is defined in Eq. (32) with the following dimensionless term:

where Vel_{surf} is found using near-surface velocity data and mathematically setting density to *ρ* = 1 kg m^{−3}, due to scaling in the dimensionless system described in Eqs. (6) and (7).

We then apply the WCE method on each wave record at truncation order K = 3, M = 2.

### A. WCE-based ocean surface, *η* kurtosis for Las Cuevas Bay

The result featured in Fig. 8 depicts ocean surface kurtosis values directly calculated at every second over a moving 60 s window (yellow line). Kurtosis is also calculated via the WCE method (blue and orange lines) with the initial condition set as Las Cuevas Bay data collected at 1350HRS or 1:50 PM, May 26, 2014, where bottom depth is 15 m and conditions are mild. Given the fact that these data are collected at 14 m depth and calculations are made via a 1D wavenumber model, ignoring directional effects, the WCE-based kurtosis estimations qualitatively resemble the directly calculated values of data kurtosis.

### B. WCE-based ocean surface energy calculations for Las Cuevas Bay

The above estimations for the *η*_{α} expansion terms from field data are used to calculate the time evolution of kinetic energy and potential energy.

In Fig. 9, we show a time series of data-generated estimations for kinetic and potential energies alongside WCE estimations for kinetic and potential energies using a 20-minute record generated at Las Cuevas Bay.

In Figs. 10–16, we show these WCE estimations for kinetic and potential energies all starting at different initial conditions defined using the various 20-minute records generated at Las Cuevas Bay. These initial conditions are indicated via spectral plots. As shown in the literature, we expect very unidirectional initial conditions to coincide with very nonlinear wave energy estimations (i.e., difficult to characterize) and very multidirectional initial conditions to coincide with more linear-like wave energy estimations (i.e., kinetic and potential energies appear to be near sinusoidal and out of phase from each other). In general, our estimates are aligned to the modeled results as shown.

### C. Discussion

The Wiener chaos expansion method allows for the calculation of the statistical moments of the ocean surface. In the present work of the wave power community, the time evolution of higher-order statistical moments may be useful to aid the estimation of the time-dependent probability distribution of the nonlinear response of a wave energy converter.^{31} Our present results address the statistical moments of the ocean surface, treating more with the environment than with the response of wave energy converter devices in the environment. These results will now be discussed.

In a linear ocean wave field, kinetic energy and potential energy are out of phase from each other. All the results for wave energy via the WCE method eventually trend toward the described linear behavior, showing a periodic exchange between kinetic energy and potential energy, although kinetic energy always appears to be less linear with higher frequency characteristics. In a realistic nonlinear ocean wave field, kinetic energy is usually greater than potential energy. This is seen in all results generated from data and from theoretical Stokes waves being propagated through the WCE system. As the numerical scheme loses energy to its numerical stabilizing process, kinetic and potential energies are numerically lost from the system and this trends the result toward linearity.

In all results, potential energy calculations in the initial condition use actual surface elevation data scaled by their standard deviation and zero-centered around their mean. As expected, the initial potential energy calculation, according to Eq. (32), is a value near 1 KJ. Similarly, kinetic energy calculations in the initial condition use actual surface elevation data scaled by their standard deviation and zero-centered around their mean, alongside surface elevation velocity data scaled by their standard deviation. As expected, the initial kinetic energy calculation, according to Eq. (35), is typically a value between 10^{0} and 10^{1}.

Throughout the eight hours of data collection, the ocean surface behavior ranged between comparatively high-energy and comparatively very low-energy, depending on the immediate environmental conditions.

Observing the results where the initial conditions have high energy, as seen in Figs. 10–12, and comparing to the results where the initial conditions have comparatively lower levels of energy, as seen in Figs. 12 and 13, and also to the results where the initial conditions have comparatively the lowest energy, as seen in Figs. 14–16, we see similar behavior in the WCE-derived values of kinetic and potential energies. In the first 40 peak periods of propagation, kinetic energy and potential energy each move through a similar number of cycles, with comparable drop-off (although less drop-off in lowest energy records) and all cases with kinetic energy showing more nonlinearity above the fundamental periodicity of the result.

We see then that the WCE-derived time evolution of ocean surface energy is stable despite differing initial conditions, confirming the suitability of the WCE method for calculating the probabilistic time evolution of ocean wave field energy. This demonstrates the possibility of the use of the Wiener chaos expansion method in understanding the probabilistic behavior of the time-evolution of ocean wave energy for ocean surface energy applications.

### D. Conclusions

Before the Wiener chaos expansion method, a universal tool suitable for deriving the probability distribution function of a nonlinear system did not exist,^{15} but now it does. The Wiener chaos expansion method is computationally inexpensive, allowing for a direct, semi-analytic calculation of the time evolution of statistical moments. In the present work, we demonstrate that the Wiener chaos expansion method is a reasonable method for calculating kurtosis and understanding the bounds of the statistical properties such as kinetic energy and potential energy of the ocean surface. Our energy results describe the energy of the ocean surface. However, our statistical moment results can be applied in the future to cumulant closure and quasi-moment closure methods in the existing literature to aid the wave energy community’s estimation of the probability distribution of wave energy converter responses to oceanic forcing.

We present our current uni-dimensional time-domain calculation that includes the effect of directionality only in the initial condition as a very strong proof of concept of the efficacy of the Wiener chaos expansion method as a deterministic means of understanding the non-Gaussian behavior of the ocean surface. We understand the further need for extending this result to a bi-dimensional time-domain calculation, and this is our intended next step. We believe that this will produce better qualitative comparisons with field-observations.

## ACKNOWLEDGMENTS

L.H. is grateful to W. Craig (deceased) of the Fields Institute, as well as L. Demanet in the MIT Mathematics Department and DKP Yue in the MIT Mechanical Engineering Department, for useful discussions in the formation of the present approach.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*et al.*, “

*et al.*, “