Numerical Considerations for Quantifying Air–Water Turbulence with Moment Field Equations

We investigate energy transfer of air–water interactions and develop a numerical method that captures its temporal variability and generates and tracks the short waves that form in the water surface as a result of the air–water turbulence. We solve a novel system of balance equations derived from the Navier–Stokes equations known as moment field equations. The main advantage of our approach is that we do not assume a priori that the stochastic random variables that quantify the turbulent energy transfer between air and water are Gaussian. We generate non-conservative multifractal measures of turbulent energy transfer using a recursive integration process and a self-affine velocity kernel. The kernel exactly satisfies the (duration limited) kinetic equation for waves as well as invariant scaling properties of the Navier–Stokes equations. This allows us to derive source terms for the moment field equations using a turbulent diffusion operator. The operator quantifies energy transfer along a space time path associated with pressure instabilities in the air–sea interface and transfers the statistical shape (or fractal dimension) of the atmosphere to the wind-sea. Because we use observational data to begin the recursive integration process, the ocean–atmosphere interaction is inherently built into the model. Numerical results from application of our methods to air–sea turbulence off the coast of New Jersey and New York indicate that our methods produce measures of turbulent energy transfer that match theory and observation, and, correspondingly, significant wave heights and average wave periods predicted by our model qualitatively match buoy data.


Introduction
In our companion paper [8], we introduce a mathematical framework to quantify airwater turbulence using recursive integration and moment field equations. The moment field equations are derived from the Navier-Stokes equations and preserve spectral moments of gravity waves of high wave number (waves with a period between 1 s and 30 s) generated by air-water turbulence. When we speak of air-water turbulence we refer to the transfer of energy between the wind and water by pressure pulsations over a duration (T < 30 s) that is shorter than the peak period of the waves that comprise the wind-sea (waves with a non-dimensional frequency ν = f m U 10 /g > 0.14). The transfer can occur from air-to-water or from water-to-air depending on the dynamics of the boundary layer.
The mathematical framework developed in [8] coalesces probability theory and observations of air-water interactions into a wave model that belongs to a class of models known as phase-averaging models; rather than attempt to quantify each individual wave train that makes up the water surface, phase-averaging models quantify the water surface in a statistical sense and predict the expected wave momentum, energy, etc., based on the physical conditions at a given geographic location. The expected wave properties are then used to predict statistical quantities of interest such as significant wave height (H s , the average of the highest one-third of all individual waves measured during a 30 min period, see [16]), average wave period, etc., which are used in wave forecasting and by engineers in design of ships and marine constructions.
The development of phase-averaged wave models is in some respects a mature field. For a concise history of our current understanding of wind-generated waves and the mathematical apparatus developed to quantify the generation and propagation of these waves, we refer the reader to [16,17]. While our collective knowledge in terms of wind-generated waves has grown over the past 60 years, the development of new algorithms that make use of this knowledge that reduce computational time in terms of computing the wave solution while still retaining accuracy has not kept pace (aside from the authors' previous work [7] and the work of Mellor et al. [29]). This is an important issue because coastal hazard management is seeking to move to an ensemblebased approach to develop resilient coastal communities in the face of an uncertain future due to climate change. Short waves that form on top of the tides (waves with a period of 12 h or 24 h) and the long waves created by the pressure drop associated with a strong storm (waves with a period greater than 2.5 min) can increase storm surge heights by up to 35% and is an important component to include when examining the risk of damage to coastal communities by storm surge [9]. However, due to the computational expense of calculating the short wave component (see [10,29]), the short waves are routinely either disregarded from ensemble simulations [24] or the convergence criteria of the solution are weakened to speed up run times [32].
Further, phase-averaged wave models typically have trouble accurately predicting peak wave period and peak wave heights in strong storms without the use of tunable parameters that modify the rate of energy transfer between the wind and water [5]. This is in part due to the fact that Mile's generation mechanism is insufficient in strong storms [5] and the calculation of nonlinear wave-wave interactions in spectral energy space broadens the spectrum too quickly (see [5,36] for further details). This results in a reduction of energy transfer from the atmosphere to the water [5].
Rather than attempt to explicitly resolve the energy transfer between the air and water in spectral space, we quantify the likelihood that energy transfer will occur between the air and water using a conditional probability that satisfies invariant properties of the Navier-Stokes equations and is a solution of the energy balance equation for the short waves. We pair this conditional probability with a recursive integration method (similar to the cascade model of Meneveau et al. [30] and the curdling model of Mandelbrot [25]) to determine how energy transfer is distributed in time. We show in [8] that using this method as input to a system of moment field equations predicts significant wave heights with similar RMS errors as compared to a third generation (phase averaging) wave model while running almost 40 times as fast. In this investigation, we show that our recursive integration methods are able to resolve peak wave heights and peak wave periods of a growing wind-sea without the use of any tunable parameters.
More specifically, what makes our approach to modeling the wind-sea unique is that the mathematical model is derived from the Navier-Stokes equations and airwater interactions are quantified using non-Gaussian measures of turbulent energy. This allows energy transfer between the air and water to occur in pulses rather than in a smooth linear fashion.
Beginning from a set of hourly averaged observational data that includes wind speed, direction, and pressure, we quantify the probability of turbulent energy transfer between the wind and water surface using a hyperbolic distribution and a recursive integration method that produces a multifractal representation of energy transfer between wind and water. The energy measure of the parameterization is intermittent, and energy dissipation on relatively small time scales (periods < 30.0 s as compared to hourly periods) is non-Gaussian (events in the tails of the distribution play a role in the process).
The multifractal representation of the energy transfer between the wind and water serve as input to the moment field equations using a turbulent diffusion operator. The diffusion operator is temporal rather than spatial, and quantifies the probability that the wave energy at a geographic coordinate, x 0 , will go from a state, n , at time t n , to a state, n+1 , at time t n+1 . The diffusion operator is a function of the speed and direction of a pulsation in the wind as well as the group velocity of the energy pulsation in the air-water interface (created by the pulsation in the wind), and the peak phase speed of the wave group that energy is being transferred to (see Sect. 2.2 for the details).
The moment field equations are a system of hyperbolic equations that generate and track short gravity wave energy on the surface of a hydrodynamic domain. The number of moments preserved in the system depends on the wave dynamics under consideration. In our initial investigation, we preserve the first two moments of the wave field which correspond to the variance of the wind-sea at a given geographic coordinate, x 0 , and the first-order perturbation from the variance, which is a measure of the characteristic wave number.
We discretize the moment field equations with a discontinuous Galerkin (DG) finite element method that puts weak constraints on the continuity of the moment field. This allows the discontinuous nature of the energy transfer between the air and water interface to be readily captured by the model. Our choice of the DG method is particularly advantageous because it allows the wind-wave model to use unstructured grids as well as local polynomial approximations of varying degree. The latter detail is particularly important when trying to resolve the peak period of a growing wind-sea. In fact, phased averaged wave models typically do not quantify the average wave period as accurately as significant wave height [16] and numerical results in Sect. 6.2.1 show that varying spatial resolution and (local) high-degree polynomial approximations are necessary to capture changes in the spectral peak wave period of a growing wind-sea.
The remainder of this paper is organized as follows: in Sect. 2, we summarize the moment field equations and discretize them in Sect. 3 with DG finite element methods. We close the resulting system of equations and outline the solution method to parameterize the atmospheric input in Sect. 4. In Sect. 5, we explore methods to generate limit distributions of the wind velocity and pressure and we present numerical results in Sect. 6. Finally, in Sect. 7, we discuss some conclusions and future work.

Mathematical Model
Let be a bounded domain in R 3 . We track energy associated with short gravity waves generated by air-water interactions within by solving a set of moment field equations, where ∇ xy = ∂/∂ xî + ∂/∂ yĵ, xy = ∂ 2 /∂ x 2 + ∂ 2 /∂ y 2 and n corresponds to the order of the moment, m n , of the one-dimensional density (F(x 0 , k, t)) of the short gravity waves at the geographic point, x 0 , In the expression above, k = ||k|| is the magnitude of the wave number vector with mean direction,θ , and k i and k f are the low and high wavenumber cut-offs of the wind-sea. The mean direction satisfies the differential equation, where D(x, t) parameterizes the inhomogeneities of the medium of propagation, see [38] for example. It can be noted that D incorporates atmospheric input, fluid depths, etc.; the c i are coefficients that represent the velocity, acceleration, etc., of the wave group (see Appendix A in [8]), and are functions of a dispersion relation, which relates the angular frequencies (ω n+1 = 2πf n+1 ) to the characteristic wave numbers,k n+1 , defined as where the bounds of integration in (5) are from k i to k f . The number of moments preserved in (1) is a function of the wave dynamics under consideration. In this numerical investigation, we set N = 1 while neglecting the effect of currents (c 0 ) and the highorder acceleration term (c 2 ), and solve the system of equations, where dE 0 /dt and dE 1 /dt are duration-limited energy pulsations created in the water surface (that travel with a peak phase speed denoted by c p f and a group speed denoted by c g f ), created by pressure pulsations in the wind; c p is the peak phase speed of the waves that the atmosphere is interacting with,Ė 0 is the energy supplied to the waves from the atmosphere, θ w is the mean direction of the wind,θ is the mean wave direction, and m 0 and m 1 represent the characteristic energy and the representative wave number of the wind-sea, respectively. Qualitatively, the term c g f / c p is similar to the inverse of the wave age χ 10 = c p /U 10 . However, the main difference between the two terms involves the fact that the term in Eq. (6) utilizes the group speed of the energy pulsation in the water surface created by a pulsation in the wind rather than the speed of the wind itself.
It can be noted that the time derivative terms (of the energy pulse supplied by the wind to the water) on the right-hand side of (6) are derived from the second-order (diffusion) term in Eq. (1) by making use of the wave equation (see [1] for example), and by applying analytical probability theory introduced by Kolmogorov [20]. More specifically, the general result we use from [20] examines the probability of a stochastic process of countable states going from a state 1 at time t = t 1 to the state 2 at time t = t 1 +dt l (where dt l is the length of duration of the transition). Kolmogorov's theorem essentially states that a n + 1-order derivative of the probability matrix of a locally continuous and homogeneous process (the process only depends on the separation in time) is equivalent to the convolution of a transition matrix and an n-order derivative of the probability matrix. In application of the theorem to our work, the transition matrix of the waves is dependent on: (1) the conservative quantity being preserved (momentum, energy, etc.) and (2) the flow of energy within the wave spectrum, which we parameterize with Hasselmann's nonlinear energy transfer mechanism [13]. The details of the derivation can be found in our companion paper [8].
It can be emphasized that the form of the moment field equations given by (6) quantifies the short wave energy of the water surface, however, the moment field equations can be used to quantify other properties of interest of the waves, such as, the wave momentum or energy flux. In this case the source term for m 0 is modified depending on the quantity being preserved. For example, in our companion paper [8], we use the energy flux form of the moment field equations to generate and propagate short gravity waves created by wind over Lake Erie.

Calculation of Mean Wave Direction
Aside from initial wind-wave generation, the direction of the wind velocity rarely coincides with the wind-sea direction, and thus, only a portion of the atmospheric turbulence will transfer energy to the short waves [16]. The rate of energy transfer is, The portion of the turbulent energy that does not actively grow the wind-sea works to turn its direction. If we integrate (7) the total amount of moment mass imparted to the wind-sea in the time interval (t f − t i ) is, which we add to the variance to determine the change in mean wave direction, i.e., δθ = tan −1 δm 0 sinθ +Ē 0 sin θ w δm 0 cosθ +Ē 0 cos θ w ,

Elucidation of the Source Terms
The source terms presented in the previous subsections quantify the energy resonance of the wind-sea; a pulsation in the energy of the wind creates a corresponding pulsation of energy in the water surface. Whether or not this pulsation grows into a longer wave depends on the speed and direction of the wind pulsation relative to the water surface coupled with the duration of the pulse. This resonance mechanism was originally posited independently by both Phillips [35] and Miles [31], however, our resonance mechanism differs in the following respects: (1) we quantify the resonance mechanism in terms of a moving frame of reference like the work of Phillips [35], however, rather than define this frame of reference in terms of a speed moving with the wind stress we define our frame of reference in terms of the wind speed associated with the flow of energy in the wave spectrum (see Sects. 4.1.1 and 4.1.2), and (2) we do not explicitly resolve the shear instability that instigates the formation of the wave in the water surface as in the work of Miles [31], but rather, we use a conditional probability to measure the likeliness that an instability (and subsequent energy transfer) will occur at the air-water interface given the current and past physical state of the wind-sea (see Sect. 4.1.1).
In fact, our parameterization of the air-water turbulence connects each state to a previous state which can be traced back in time to a generating event. In the case of airwater turbulence, the generator of the wind turbulence could be due to the wind flowing over rough topography and a rough coastline before interacting with the water surface. The distribution of eddies (or regions of energy dissipation) created by interactions with these rough boundaries does not occur everywhere within a fluid but only occurs within a fraction of it. The dimension of the space where this energy transfer occurs is known as the fractal dimension. Each interaction injects a different statistical pattern (or fractal dimension) of energy transfer into the fluid flow. The result is that the fluid flow contains a spectrum of fractal dimensions of energy transfer and the flow is said to be multifractal [12].
Because we are interested in small-scale turbulence (turbulence with a period smaller than the peak period of the waves), we utilize a recursive integration method (see Sect. 5.2) to create a refined temporal distribution of the wind (from hourly averaged observational data) that is multifractal. The recursive integration proceeds in a similar fashion to other cascade models such as the p model [30], β model ( [12,25]), and multifractal model ( [12,25]), however, rather than use a priori considerations (in terms of the refinement weights) to refine the spatial distribution of the energy transfer our recursive integration methods make use of the current and past states of the flow field to determine how to refine the temporal distribution of the kinetic energy (see Sect. 4.2). This is made possible by the fact that for fluid flows with small viscosity the Navier-Stokes equations are invariant under scaling transformations in the form of self-affine power laws (see Sect. 4 and [12]). We leverage this fact and model the wind and water as an inviscid fluid that contains singularities (which is where the initial energy transfer occurs), and we incorporate viscous effects (dissipation) via a closure scheme that has been verified by direct numerical simulations (DNS). More specifically, the exact quantitative relationship between the probability of energy transfer and subsequent growth of the wave is based on the original energy transfer mechanism of Hasselmann [13] which takes into account atmospheric input, nonlinear wave interactions, and dissipation and was verified by the DNS of Tanaka [39] (see Sect. 4.1.2).

Numerical Discretization of Moment Field Equations
Our mathematical model consists of a system of balance equations that preserve the spectral moments of wave properties (momentum, energy, or energy flux) of the windsea. Rather than discretize the strong form of these equations, we recast them into a weak form that admits discontinuous solutions-a necessary measure due to the intermittent nature of air-water interactions-see Fig. 3 for example.

Weak Forms and Subspaces
Given a two-dimensional partition of the air-water interface (of the hydrodynamic domain ) into a set of non-overlapping elements (or cells) e denoted by, T h , we obtain a weak form of the moment field equations if we first multiply (6) by a suitable smooth test function υ ∈ V, integrate over each element e ∈ T h , and integrate the flux term by parts to obtain, whereê is the outward unit normal of the element boundary ∂ e and n = 0, 1. Setting υ = 1 in (9) gives the standard weak form of the monotone finite volume method where a large amount of literature exists in terms of numerical solution methods, see [23] for example. In this initial investigation we utilize a discontinuous Galerkin (DG) method [6], and set the finite-dimensional space of test and trial functions to where P p demarcates the space of polynomials of at most degree p. Now, rather than solve (9) directly, we search for trial solutions from the finite-dimensional space of test functions defined by (11). Specifically, given a set of basis functions φ = [φ 0 , φ 1 , . . . , φ κ ] T , we express the trial solutions as (1) , θ (2) , . . . , θ (κ) are the time dependent degrees of freedom of the finite element solution. It can be noted that we use products of Jacobi polynomials of degree κ, P p κ p=0 , as our basis for V hp . The orthogonal triangular basis is defined in terms of a "collapsed coordinate" system and results in a matrix free implementation of the method, see [21] for more details.

Numerical Flux and Discrete Weak Form
The finite element space defined by (11) does not explicitly enforce continuity across element boundaries, and therefore, to ensure a consistent flow of information throughout the spatial domain we replace the dual-valued flux in (9) with a so-called numerical flux that is a function of left and right limits of the trial solution. For example, given an arbitrary function υ h ∈ V hp at an element boundary point, x i , we set the left and right limits of the function to respectively. In this particular work we replace the dual-valued boundary flux with a local Lax-Friedrichs (LLF) flux, namely, is the maximum peak phase speed of the short waves and n = 0, 1. The discrete weak form of the problem is then: find m n h ∈ V hp such that for all test functions υ h ∈ V hp the expression, holds over each e ∈ T h with n = 0, 1. Application of the DG spatial operator results in a system of ODEs where L h corresponds to the DG spatial operator andL h is a function of the inhomogeneities of the medium which reduce to simple vector additions in the case of the wind-sea, see Sect. 2.1. We evaluate the integrals in Eq. (12) using numerical integration rules of sufficiently high degree, and discretize (13) with so-called strongstability-preserving (SSP) Runge-Kutta (RK) methods [22], where for each time step (from t k to t k+1 ) we have, 3. Finally, set U It can be noted that U h = [m n , θ] T and h is a slope limiter that stabilizes the method when sharp gradients form in the numerical solution for p > 0, δ j t is a sub-time step of the time step t, and the α i j and β i j are coefficients that define the order of the time stepper, see [22] for more details.

Numerical Parameterizations
To close the discrete DG system of Eq. (12), we need to determine the characteristic constants that define the air-water turbulence over a given temporal region denoted by R j (which is defined in the next section).

Quantifying Atmospheric Velocity and Pressure
To quantify the temporal distribution of small-scale wind turbulence that generates the wind-sea, we consider the Navier-Stokes equations, where p is the pressure, ρ a is the density of air, and μ is the kinematic viscosity. Frisch [12] has shown that in the limit μ → 0 the Navier-Stokes equations are invariant under scaling transformations The scaling also holds for small viscosity if [12], This implies that the dissipation (or energy transfer), defined as = μ (∇u) 2 ( · denotes ensemble averages), scales as [12] → λ 3h−1 .
From this it follows that the absolute probability of encountering an active region of turbulence is, for η < r << L, where η is an inner cut-off scale and L is the integral length scale, see [3,12] for more details.
Kolmogorov's seminal 1941 paper [19] (one of the most influential papers in terms of scaling laws and turbulence) assumes that nonlinear interactions are local in momentum space (k-space). This implies that energy dissipation occurs everywhere within the fluid and corresponds to a universal value of h = 1/3. However, experimental data and observations suggest that h is not a universal constant (particularly for geophysical flows), but varies within a given fluid domain ( [4,12]). A number of theories account for this fluctuation in h, one being the notion of multifractal turbulence ( [12,25]). However, this theory, for all of its merits is mostly phenomenological. It generally lacks a sufficient way to generate the distributions of h aside from a priori considerations. The interpolation method we present in Sect. 5, while multifractal, naturally generates the distributions of h by preserving large-scale averages and end point relations defined over a set of temporal local regions R n j j=1,M n=0,N that pave a given space-time curve fixed at the geographic point x j . Application of the method to observational wind data collected by the National Data Buoy Center ( [33,34]) results in scaling exponents h ∈ [1/3, 0.42], see Sect. 5 and [8] for the details.

Conditional Probability and Self-Affine Velocity Kernel
Rather than utilize the absolute probability (17) that is highly dependent on L, we utilize a conditional probability to generate measures of the velocity and energy transfer connected to air-water interactions [8]. That is, given the condition that an active region of turbulence occurs within the period l d η << l d < t l < l 0 << L T (where l d η is the molecular length scale, l d ≈ 1 × 10 −6 and L T is the temporal integral length scale), we assume that the local distribution of the atmospheric velocity at the spatial point x 0 is temporally hyperbolic, with magnitude [8], and direction, where t l is the local time, g is acceleration due to gravity, (||u 0 || , θ 0 ) are local constants, and ω θ is the minimum turning frequency of the characteristic reference velocity (u 0 ), i.e., the non-pulsating component of the large-scale flow field of the space-time region under consideration. It can be noted that the exponents (q, ϒ) quantify the velocity pulsation, u, over the time interval l 0 − l d from the reference velocity field u 0 (if q = 0 and ϒ = 0 then u = u 0 ). The union of all the local regions used to pave the temporal velocity curve gives the global distribution of the wind velocity. For a given time series we have, where the magnitude ( u ( j) ) and direction (θ w ) are functions of the spatial coordinate, x 0 , and global time t. It can be noted that N is the total number of local regions used to pave the temporal wind velocity curve and the l attached to u l indicates that the continuity of the temporal representation of the distribution of the wind velocity depends on the inner cut-off of the local distribution. In particular, we define a temporal region R ( j) located at the geographical coordinate x 0 with height, Each R j is a local reference frame for the dynamical region under consideration, and is characterized by a set of characteristic constants u Each R ( j) is a local reference frame (moving with the constant speed ||u 0 ||) that contains the dynamic fluctuations of the flow field in the time interval l 0 (at the point x 0 ). Embedded within each R ( j) is a cut-out region ]l d , l d [ where functions (18) and (19) are ill defined, see Fig. 1. This region corresponds to a measure of uncertainty in terms of the measure of the turbulent velocity on the length scale l 0 (a singularity occurs within this region). Numerically, the requirement for a consistent (and stable) integration of (13) is that l d < dt (where dt is the time step of the numerical integration). This ensures that the turbulent energy transfer occurs over a period that is less than the period of the waves. We examine this issue in more detail in Sect. 5. The total pressure distribution ( p ( j) T = p 0 − p ) associated with (4.1.1) over a given R ( j) written in terms of the average pressure,p, over a given duration, is derived in [8] as, where the global distribution of the total pressure is the union of the pressure kernels associated with each R ( j) , see [8] for the details.

Duration-Limited Solution to Wave Energy Balance
To close the moment field equations, we need to quantify the energy resonance of the wind-sea. In this work, we leverage the fact that closed form solutions exist to a parameterized form of the energy balance equation for waves when the magnitude of the wind can be represented locally by (18). In other words, by representing the wind turbulence as a series of pulses (given by (18)), we immediately have access to mathematical relations that quantify the corresponding pulse of energy that forms in the water surface. More specifically, Hasselmann [13] derived a rigorous parametrization of wind-sea growth based on observational and experimental data that takes into account wave growth via atmospheric input, nonlinear wave-wave interactions, and dissipation. The energy transfer mechanism and the corresponding flow of energy in the wave spectrum has been verified by the DNS of Tanaka [39]. The parameterized energy balance equation Hasselmann developed typically must be integrated numerically, however, when short gravity waves are duration-limited (as the turbulent pulsations in our model are), then the non-dimensional representative frequency (ν = f m U 10 /g) and non-dimensional energy ( = E g 2 /U 4 10 ) of the corresponding wind-sea can be expressed solely in terms of the characteristic constants (see [13] for details), where and The shape parameter, Λ, relates the non-dimensional energy, frequency and dissipation parameter, α ( j) , via where Λ = 1.60(±0.02)(10 −4 ) arises from wind-sea observations [13]. The measure of the high-frequency cut-off, α ( j) , is expressed in terms of u ( j) 0 and q ( j) as, where b ( j) is given by (24), and ν ( j) by (23). It can be noted that the high-frequency cut-off is related to the high wave number cut-off via the dispersion relation (4). The global distributions of the non-dimensional frequency and energy of the durationlimited wind-sea are, Because the wave pulsations that form in the water surface are of a much shorter wave length than the depth of the water, the magnitude of the group velocity, ||c g f || , of the wave pulsation is, and the peak phase velocity is, We use (23) to determine the dE 1 /dt terms in the moment field equations. Making use of (23) and differentiating with respect to time yields, where , and C 0 = 5.25 10 −4 and C 1 = 4.20 10 −3 . It can be noted that in some instances (e.g., in the absence of appreciable advection) Eq. (23) can be used solely with our recursive integration schemes (to be defined in the next section) to capture the variance of the wind-sea to predict significant wave height. In this particular case, the variance of the water surface can be calculated via (27) and we define the significant wave height as, where H is the Hurst exponent and is equal to H = 2 − D G , where D G is the graph dimension calculated via the variation method. (In practice, an approximation of H can be obtained by setting H = 1/D where D is the dimension of the support of the energy transfer, see Sect. 5 and [8].) In the wave literature, H is set equal to 1/2 because the distribution of small amplitude short waves is close to Gaussian, see [16] for instance.

Determination of the Characteristic Constants
We determine the set of constants {||u 0 || , θ 0 , q, ϒ} that characterize the local distribution of the wind velocity through preservation of large-scale averages combined with some reconstruction of the end point values of a given R ( j) , see Fig. 2. More specifically, for a given R ( j) fixed at given geographic coordinate, x 0 , of length l 0 the initial and final values of the magnitude of the velocity are, Solving for q ( j) we have, Setting these equations equal to one another and simplifying leads to an expression that can be solved for u The exponent q ( j) is determined through the preservation of large-scale averages of observational wind data, i.e., the area. In particular, the integral of (18) is where t l is the local time,ū ( j) is the average wind speed over the region R ( j) . We solve (34) and (35) simultaneously with the help of a Newton solver [1]. More specifically, given a system of equations of the form with an initial iterate in the neighborhood of the true solutionx, i.e.,x = x (0) + x, we can expand (36) about x (0) using a multidimensional Taylor series, and reduce the nonlinear system to a linear system that can be solved for x [1]. In matrix form the solution iteration becomes where J −1 is the inverse of the Jacobian matrix, Because the velocity kernel is analytic we can evaluate the Jacobian and its inverse exactly, where dt = l 0 −l d . In some instances u i or u f can take on a value close to zero causing the Jacobian matrix to become singular. In order to circumvent this issue we utilize a shift constant on the end point relations of the velocity, i.e., we set u i = u i + C s and u f = u f + C s where C s = 1.

Limit Distributions of Atmospheric Input
We now apply the nonlinear solution method outlined in the previous section to observational wind data obtained from NDBC buoy 45005 in Lake Erie [34] and examine the nuances of canonical versus microcanonical recursive refinement. The goal is to produce a fractal representation of the wind velocity whose dimension matches the observations of Syu et. al. [39] and Richardson [37] (as analyzed by Hentschel and Procaccia [14]), which place the graph dimension of the atmospheric turbulence at D G = 1.60 ± 0.03.
The graph dimension is directly related to the self-similar fractal dimension, D, in the unifractal case [26]. More specifically we have, so that a graph dimension value of 1.60 gives a self-similar fractal dimension of D = 2.50. In the expression above, H is the Hurst exponent, and is important in a statistical context because it gives a measure of the memory of an observable process whose average is zero. For instance, say we want to observe the movement of a particle in time (denoted by B H (t) [26]) and the generator of the particle's motion is unknown to us. If we observe that the average movement of the particle's motion measured in reference to a zero datum is null, and the energy is finite, then the expected variance of the particle's motion over a duration of time T can be expressed in terms of the Hurst exponent, [26], where V[·] is the expected variance and 0 < H < 1 [26]. Setting t = 0, the correlation between a past duration (−B H (−T )) and a future duration (B H (T )) is also expressed in terms of H [26], Dividing the above result by the variance (T 2H ) gives the correlation [26], 2 2H −1 − 1.
If H > 1/2 the correlation is positive and the motion of the particle is persistent; a given trend in the motion of the particle will persist for long periods of time [26]. If H < 1/2, then the correlation is negative and the motion of the particle is said to be anti-persistent; the motion of the particle has a tendency to return to some past state in time [26]. In the case that H = 1/2 the motion is "memoryless", and the variance is connected to the standard deviation via the square root operator [26]. It is in this case, and this case only, that the motion truly can be deemed as random; in the case of either persistent or anti-persistent motion the movement of the particle is not random, but rather, the generator of the motion is opaque to the observer. It can be noted that the Hurst exponent (H ) is related to the scaling exponent (h) of the Navier Stokes equation via, For example, Kolmogorov's 1941 theory of turbulence assumes h = 1/3. The Hurst exponent equals 1/3 and the motion of the turbulent fluctuations is anti-persistent. The self-similar dimension of the energy transfer is D = 1/(1/3) = 3 and occurs nearly everywhere within the fluid, while the graph dimension of the turbulence is D G = 1.67. As stated previously, however, atmospheric observations tend to show that the graph dimension of the wind is less than 1.67 [39]. The turbulence is intermittent and does not occur everywhere within the fluid. Further, H is not a constant but varies over a given range of values where each value corresponds to a different generator of turbulent energy transfer within the fluid, see Tables 3, 4, 5, 6 where H = 2 − D G , for instance. The fact that H < 1/2 for turbulent flow is a result of the vortexes created by the generator of the turbulence.
To study the graph dimension of our parameterization of the wind velocity, we first use the methods presented in Sect. 4.2 to parameterize the wind velocity with a local reference frame (the R ( j) ) that are the length (or duration) of the hourly averaged wind data, i.e., l 0 = 3600 s. We then calculate the resulting graph dimension of the velocity distribution of the wind using the variation method [11].

The Variation Method
We follow the discussion in Dubuc et al. [11] and Syu et al. [39]. Let f (t) be a graph, and let it be covered with N equally spaced regions of length, The total variation of f (t) is defined in [11] as, which is the total area of the regions required to cover f (t). In other words, V is a measure of the variation of the space-time region under consideration-the rate of the In practice, D G can be calculated through a least squares approximation of the points as ε → 0, see [11,39] for more details.
Using the variation method we calculate the inner cut-off of the wind velocity necessary to produce a graph dimension D G ≈ 1.60 for NDBC buoy 45005 [34] for the years 2010 and 2015. It can be noted that the inner cut-off in this case corresponds to a relative zero for the variational method. The results are shown in Tables 1-2 where the salient points to make note of are: i. in both 2010 and 2015 the inner cut-off of the turbulence appears to be a seasonal phenomenon and ii. the inner cut-off of the turbulence (for hourly averaged data) ranges from 56 s down to 17 s, values that can be too large to consistently integrate the moment field equations with explicit SSP Runge-Kutta time steppers (due to the fact that the time step is limited by a CFL condition).
We remedy this issue by using a recursive integration process to shrink the region of uncertainty associated with each velocity distribution u ( j) in a manner such that ]l d , l d [ << dt. We proceed by first subdividing each R ( j) into two subregions R ( j) 1 and R ( j) 2 ; we then calculate the expected velocity of each subregion and pair these values with the nonlinear solution method described in Sect. 4.2 to create a new set of velocity kernels. This new set is composed of twice as many R ( j) as the original set and adds cut-outs along the space time path while shrinking the region of uncertainty ]l d , l d i [ associated with these cut-outs toward zero.

Recursive Integration and Velocity Measure
Consider the discrete sequence of lengths associated with a given R ( j) , The average, or expected speed of each subregion of level i + 1 is connected to the expectation of level i viaū where C i+1 is the number of subregions of level i + 1 and the w ( j) s correspond to the measure (or weights) of each subregion s, which are defined as, The expectation of the parent kernelū The offspring,ū ( j) s , corresponds to the expected speed of the turbulent velocity over an interior time interval, When dyadic intervals serve as the refinement mechanism, the averages used to initiate the next level of integration are, It can be noted that the interval ]l d , l d i [ is cut-out of the integration. The manner in which this interval is handled during the recursive integration leads to markedly different measures and dimensions of the support.

Conservative Recursion
In the microcanonical approach, we set the cut-out lengths to ]l d , l d i [ = 1 × 10 −7 and hold them constant with each level of integration. The sum of the weights in (45) equals C i+1 and the refinement cascade is said to be conservative [25]. In this case ]l d , l d i [ is always less than dt, however, the graph dimension of the resulting velocity distribution will not be in a range that matches theory and observation until after five levels of recursive integration. The dimension of the support in this case is approximately equal to the embedding dimension E (the resulting measure of the turbulence is E-filling). In this case, the recursive integration does not dissipate any velocity to the background flow. The resulting graph dimension is D G = 1.64 ± 0.03 (D = 2.78), a similar value to many studies of fully developed turbulence (D = 2.75) as well as Kolmogorov's 1941 theory (D G = 1.67).

Non-conservative Recursion
In the non-conservative approach, we allow the cut-out ]l d , l d i [ to vary in length (more specifically to shrink or become thinner) with each integration. In this case the sum of the weights only equals C i+1 on average, local conservation is not maintained and the recursion is said to be canonical [25]. In this scenario, the graph dimension of the velocity distribution is in the range of observations and theory with each recursion and we can re-write (44) as is the value of velocity cut-out from the refinement of level i. Numerical experiments reveal that once the initial inner cut-off is established for the longest R ( j) (in this case R ( j) = 3600 secs), then subsequent inner cut-offs for each recursion follow the relation l d i+1 ≈ l d i /14. We choose the initial inner cut-off so that the normalized structure functions of the velocity distribution of the boundary layer turbulence match theory and experiment, see [8] for more details. The resulting graph dimension of the wind velocity calculated for different months in the years 2010 and 2015 are shown in Tables 5 and 6 where the values range from D G = 1.57 to D G = 1.63.

Variation Past the Inner Cut-off
If one examines the variation of the distribution of the boundary layer turbulence inside the local inner cut-off l d it vanishes. This was predicted by Mandelbrot [25] and is reflected in our numerical experiments. Once the length of the sampling regions (2 ) of the variation method are smaller than the length scale of the inner cut-off, the fractal dimension goes to a value of 1 (as → 10 −7 ), or in other words, the variation about the fixed points goes to zero and the energy associated with the fluctuations becomes null.

Numerical Results
We evaluate our recursive integration methods and the moment field equations in terms of generating and propagating short gravity waves in the water surface. In particular, we apply our recursive integration method to data collected from NDBC Buoy 44065 [33] located off the coasts of Long Island and New Jersey and model significant wave height during the later days of October 2012 during the post-tropical cyclone event known as Superstorm Sandy [2]. We also verify our DG numerical solution method against an analytic test problem put forth by Hasselmann et. al. [13], and examine the numerical recipe necessary to properly resolve the peak wave period of a growing wind-sea.

Superstorm Sandy Significant Wave Height
At the end of October 2012, a large post-tropical cyclone known as Superstorm Sandy made landfall in parts of New Jersey and New York. Storm surge reached as high as 2.74 meters in Staten Island and Manhattan and caused an estimated 70 billion US dollars in damages [2]. Because the storm was one of the costliest in U.S. history, we apply our recursive integration methods to observational data obtained from NDBC Buoy 44065 for the dates corresponding to the storm, and hindcast significant wave height using relation (31). Figure 3 shows the distribution of the magnitude of the wind velocity (the red line) generated from the hourly averaged data recorded at Buoy 44065 (represented by the blue dots). It can be noted that 4 non-conservative recursions were used to generate the velocity distribution. The dimension of the self-similar spatial trail The max and min q values are 20 and 13 times greater than the standard deviation of q (0.0017), signaling that they are statistically rare events, see Fig. 4. Finally, Fig.  5 displays modeled significant wave height using our recursive integration method coupled with relation (31) versus observed significant wave height. It can be noted that the numerical results are averaged over a 60-min window and match the NDBC buoy observations qualitatively well when wave advection is negligible in the short wave field (which occurs approximately from the 25th of October to the 29th of October). Figure 6 shows numerical results for the average wave period generated from our recursive integration methods using (23) and setting the average wave period (T wave ) equal to,T Results in Fig. 6 also match observations qualitatively well within the days 25 to 29 th when the waves are mostly wind generated. However, outside this period, numerical results for mean wave period deviate more from observations than for significant wave height, indicating the presence of combined sea and swell, and demonstrating the need to take into account wave advection in order to properly resolve the short wave dynamics.

Hasselmann's Idealized Wind-Sea
We verify our DG numerical solution method of the moment field equations against the steady-state solutions of Hasselmann et al. [13]. The test case consists of an idealized wind field blowing over a large open body of water. The wind field varies spatially according to the power law, where χ is the fetch coordinate, χ ∈ [0, L χ ], L χ is the fetch length, u χ is a constant wind speed, and p is a growth/decay factor of the wind speed over the fetch. Exact solutions exist for the non-dimensional representative frequency, variance, and dissipation cut-off of the short gravity waves (see [13]), = Cν −10/3 , At χ = 0 the solutions satisfy an initial condition of vanishing energy, i.e., = E ||u|| 4 /g 2 = 0. Note that the wind field is a homogeneous function of χ , i.e., it only depends on the separation χ − χ . We utilize the methods of Sect. 4.2 to construct a wind speed representation that corresponds to (47) spatially, where we initially set u χ = 2.0 m/s and p = 0.10 in Eq. (47). This corresponds to a nonuniform wind velocity that grows with fetch length, see Fig. 7.
Results are shown in Figs. 8, 9, where we set m 0 = E and m 1 = f m m 0 . It can be noted that using P 4 polynomials with 2 elements (10 degrees of freedom (dof)) gives a lower L 2 -error measure (||m 0 − m 0h || 2 = 7.4 10 −4 ) than using P 1 polynomials with 16 elements (32 dof, ||m 0 − m 0h || 2 = 0.0022), which demonstrates a potential benefit of using high-order local approximations in the wind-sea (lower errors via fewer degrees of freedom). It can be noted that the L 2 -error measure is defined as, where m n is the known solution and m n h is the numerical solution. Finally, regardless of the polynomial approximation utilized in (12) the numerical solution is first-order accurate, and is most likely due to incongruities in the wind interpolation. An important nuance of this particular test case worth drawing attention to involves the fact that we must time march the numerical solution from an initial condition of zero energy at t = 0 to the steady-state solution that occurs at t = T s . More specifically, when t = T s , then q = 0 in Eq. (18) and ||u|| = ||u 0 (x)||, and the moment field source terms no longer depend on time. This is implemented numerically by setting t l in Eqs. (7) and (8) to t l = T s for t ≥ T s , where T s is defined so that Eq. (23) equals Eq. (49), see Appendix A for the details.

Frequency Calculations
Historically, short gravity wave models have had trouble calculating the representative frequency (or period) of the sea-surface [9,18]. A number of conjectures exist as to why this problem persists and ideas range from a misrepresentation of wind-sea-swell interactions to error in atmospheric input. In hopes of shedding light on this issue, we We begin by setting u χ = 10 m/s and p = 0.05, which once again corresponds to a wind-sea that grows with fetch. Examination of the results shown in Fig. 10 immediately reveals that even though we capture the zeroth order moment well, we have trouble capturing the frequency in the spatial region where the wind-sea is rapidly growing. In fact, when we use P 5 polynomials with 16 elements we still significantly miss the frequency, f m = m 1 /m 0 , ( f m − f m h 2 = 1.955) even though we capture the variance well ( m 0 − m 0 h 2 = 3.93 10 −5 ). In an effort to remedy this situation we switch from a uniform distribution of elements to a geometric distribution, where x j = γ j x 0 , and try to resolve the spatial region where the peak wave frequency is rapidly changing. We initially set γ = 83/80 and x 0 = 50 and utilize P 1 polynomials to calculate solutions for the variance and frequency. As seen in Fig. 11 the numerical approximation significantly overshoots the frequency near χ = 0 even though we use a minimum element size of 50 meters and 63 elements to discretize the domain. In fact, it is very difficult to find values for x 0 and γ that allow the frequency to be accurately captured with P 1 polynomials near χ = 0.
The situation becomes even worse if we increase the value of p in (47) (the L 2 -error measure in terms of frequency increases from a value of 2.133 to a value of 4.208 when p = 0.10). This is due to the fact that the exponent p controls the growth of the frequency near χ = 0, i.e., the larger the value of p the larger the jump in the frequency. However, we can remedy this situation if we couple the geometric distribution with a high-degree polynomial basis. In particular, if we set γ = 4, x 0 = 50, and utilize P 5 polynomials we can accurately resolve the sharp change in the frequency near χ = 0 with only 3 elements, see Fig. 12, which demonstrates the power of using a nonuniform mesh with high-order polynomials. These test results indicate that even if the error in the atmospheric input is low, large errors can still accumulate in the frequency calculation if adequate resolution is not utilized in regions where the atmospheric velocity is rapidly changing. This applies in terms of both the spatial discretization and numerical discretization. These results are in accordance with Cavaleri [5] who has shown that phase-averaged wave models typically underpredict peak wave heights and peak wave periods in strong storms unless "strong, but effective tuning" is used.
Another point worth drawing attention to is the connection between the p parameter in Eq. (47) and the q parameter in the velocity kernel (18). In particular, both measure the strength of the energy transfer between the water and atmosphere. (In the case of the wind-sea test case the maximum q value when p = 0.05 is q = 0.1823.) Knowing that the q parameter is an indicator of the strength of the energy transfer, we can perhaps leverage this information to better resolve large changes in the wind velocity (which is a function of q and ϒ). For instance, when either of the exponents q or ϒ in Eqs. (18) and (19) take on values greater than a prescribed value then we can use adaptive mesh refinement (h refinement) and polynomial enrichment ( p refinement) to dynamically resolve these spatial regions where the wind velocity is rapidly changing. This could be implemented in practice by setting a threshold for triggering h and/or p refinement in terms of the standard deviation of q or ϒ.

Conclusions
We have detailed the numerical considerations necessary to quantify air-water turbulence with moment field equations and recursive integration. Our methods decompose air-water turbulence into a series of pulsations and quantify the probability that energy transfer will occur between the wind and water surface (during the pulsation) via a conditional probability. More specifically, given the condition that energy transfer occurs during the pulsation, the magnitude of the energy transfer is measured by the power law exponent q. The larger the absolute value of q the greater the probability for the transfer of energy between the wind and water surface during the turbulent wind pulsation.
It is easy to identify these regions of energy transfer in the graphs of q and the wind speed as shown in Figs. 3 and 4, where the shape of the pulse (or the distribution) of the velocity over a given local duration, l d < t l < l 0 , can be thought of as follows; if someone samples the value of the wind velocity over the period, l d < t l < l 0 , most of the change in the velocity will occur over a brief period of time and remain relatively steady before and after this change (the change in the velocity in time will not be constant as in a linear approximation). Further, the sampled velocity values will be near the average value of the velocity of the duration and the distribution of the velocity values about this average will be nearly Gaussian (especially when q is close to a value of zero), however, it should be emphasized that near Gaussian is not synonymous with Gaussian (see Frisch [12] for a discussion of this point with respect to turbulence and Table [40] for some statistical consequences of using Gaussian approximations for near Gaussian phenomena). In the case of air-water interactions, linear interpolation of the wind speed will not properly discretize the geometry of the energy transfer between the wind and water especially when interpolating to smaller scales where the energy transfer becomes intermittent. Of course linear interpolation can be used along with an appropriate "fix" (usually some tunable parameter that is problem dependent) that adjusts the energy transfer between wind and water, however, we believe it should be the goal of the mathematical sciences to move away from these "fixes" whenever possible. (This, in large part, is why machine learning is so valuable because it allows for nonlinear interpolation of data.) Interestingly, observations of water surface heights indicate that the water surface height distribution is near Gaussian, just like the wind velocity, and is a reason why the typical stochastic methods (see [16], for instance) used to fully discretize spectral space in third generation wave models have been successful. However, because of some of the simplifying assumptions necessary to make the computations more tractable, there can be a misrepresentation of the flow of energy within the wave energy spectrum (for example, the discrete interaction approximation (DIA) typically used to quantify nonlinear wave interactions forces the wave spectrum away from its true equilibrium form [36], while the use of the phase-amplitude model and the Gaussian assumption of stationarity imply that individual waves do not interact with other waves [16]). Resio et al. bring attention to this issue in [36] and they develop a list of tests that phase-averaging wave models should meet in terms of energy transfer. This is in part why we choose to use the k-space version of Hasselmann's spectral energy transfer scheme [13] because it has been verified by direct numerical simulation [39] and tightly couples the moment field source terms to the wind velocity.
One can think of the moment field source terms in the following fashion; rather than attempt to model and discretize the flow of energy in the full wave spectrum using simplifying approximations, we quantify the flow of energy in terms of the moments of the energy density via a nonlinear parameterization. The parameterization calculates the probability of encountering a region of turbulence (and energy transfer) in the air-water surface via a conditional probability that satisfies invariant properties of the Navier-Stokes equation and the nonlinear energy balance equation of the short gravity waves. Because we assume that the energy transfer occurs via a resonance mechanism, the source terms are functions of the exponent, q, as well as the parameterized wind speed, ||u 0 ||. The source terms quantify the rate of change of energy within the wave spectrum due to atmospheric input, nonlinear wave-wave interactions and dissipation and are based on the original energy transfer mechanism of Hasselmann which has been verified by the DNS of Tanaka [39]. It can be noted, however, that the source terms developed in [8] and this work are not synonymous with the source terms of Hasselmann presented in [13]. While our formulation does collapse down to Hasselmann's in the purely (steady) fetch-limited case, the source input into the moment field equations is inherently different due to the way that we decompose energy transfer between the wind and water into turbulent pulsations, as well as due to the fact that we include a pseudo-wave age term in the formulation.
Further, while we have made use of Hasselmann's energy parameterization [13] for the model closure, the mathematical structure of the moment field equations is general enough that alternative closure schemes can be used, such as, for instance, the Zakharov-Resio-Pushkarev (ZRP) wind input source term (Zakharov et al. [40]). The key difference between the two closure schemes involves the fact that the ZRP terms do not include the parameter, q. This means that the characteristic wind speed, u 0 , simply corresponds to the average wind speed over the observational period (the wind speed over a given pulsation is steady (q is a measure of the unsteadiness over a given pulsation)), and thus, to refine the distribution of the wind speed to smaller scales the recursive integration algorithm in Sect. 5.2 would need to be slightly modified. In this case an additional measure would need to be introduced to take the place of q, say, a function that measures the difference between the current wind speed and the previous wind speed (which to some extent is what q measures). A drawback of the ZRP closure scheme lies in the fact that it requires a nonlinear solve (i.e., a Newton iteration) to balance the flow of energy in the wave spectrum, however, the scheme is more general in the sense that it can be used in a phase-averaging model that explicitly models the flow of energy in the wave spectrum via the action balance equation. In this case, one can use the recursive integration method presented here to interpolate U 10 to smaller time scales and couple this with the ZRP input source term to capture the non-Gaussian behavior of the energy transfer in a fully spectral, third generation wave model. This would proceed in the following fashion; first U 10 would be interpolated to small time scales using 4-7 cycles of the recursive integration, and then this result would be averaged over a time window corresponding to the time step of the numerical integrator of the action balance equation, where the average values associated with the time step dt corresponds to the expected value of the wind velocity on this time scale.
Results from Sect. 6.1 indicate that the recursive integration methods detailed in Sect. 5.2 produce significant wave heights that match observational data qualitatively well when wave advection is negligible. When the full system of moment field equations must be solved (because wave advection is non-negligible), an attractive feature of the structure of the equations is that they are hyperbolic, and are amenable to advanced numerical discretization techniques such as the discontinuous Galerkin (DG) method (which allows for the use of locally high-order polynomial approximations) and adaptive mesh refinement. Results from Sect. 6 indicate that both high-order (local) polynomial approximations and varying spatial resolution are necessary to accurately quantify the wave period of a growing wind-sea.
Future work will focus on integrating adaptive p and h refinement as well as incorporating wind-sea-swell interactions, swell dissipation, and coastal effects into the moment field equations. Swell dissipation, in particular, remains a large source of error in phase-averaging wave models (this is in part due to a lack of physical understanding of the mechanism that leads to swell dissipation), and we plan to work with physical oceanographers to better quantify this process as well as incorporate coastal effects potentially using source terms similar to Holthuijsen et al. [15].

Data availability
The data used in this investigation can be obtained at https://www.ndbc.noaa.gov.

Conflict of interest The authors report no conflict of interest.
Computer code availability The computer code used in this investigation can be obtained at https://github. com/coltonjconroy.

Appendix A: Steady Energy Matching Condition for Duration-and Fetch-Limited Waves
In the case of a purely steady wind (q ≡ 0) blowing over a fetch-limited body of water, we ensure the duration-limited source terms supply a consistent amount of energy to the moment field by setting Eq. (23) equal to (48), and solve for t l , which yields, We then set t = t l in the source terms in Eqs. (7) and (8). There are two important points to make note of i. t l is a function of the fetch coordinate, χ ∈ [0, L χ ], where L χ is the total fetch length, and ii. in Hasselmann's parameterization, the energy of the water waves is a function of the wave frequency, and this is why we choose to match the frequency of the duration-and fetch-limited solutions (this also conveniently ensures that the first-order moments will match as well). It can be noted that in simulations of physical wave records, q is rarely (if ever) exactly equal to 0. In fact, in our companion paper [8] we model surface wave heights over the fetch-limited domain of Lake Erie and q never equals a true zero and model results match observations quantitatively well.