The procedure for deriving reduced equations by expressing the dependence of the dependent variables on a certain, geometrically distinguished spatial coordinate analytically is widespread in nonlinear pattern formation and wave dynamics. As a first attempt, one can consider the Shallow-Water or Saint-Venant equations  , which are systematically derived from the hydrodynamic basic equations for an inviscid potential flow by solving iteratively the Laplace equation for the velocity potential function. Thereby, all expansions are controlled by a geometry parameter
where d is the layer’s depth and a typical lateral length scale (e.g. wave length). A similar expansion can be done for highly viscous problems and the so-called Thin-Film equation is obtained which describes the dynamical behavior of the shape of the film’s surface    . For convection problems, a projection onto certain z-dependent modes allows for the reduction to 2D systems where only the horizontal coordinates occur. Standard equations in normalized form like the Allen-Cahn, the (complex) Ginzburg-Landau or the Swift-Hohenberg equations can be derived, partially by heuristic arguments concerning at least their nonlinear parts  .
Stokes  initiated the theory of nonlinear deep-water waves. Three branches of development for strongly nonlinear theory have grown from seminal papers. 1) The highest Stokes wave  . 2) The instability of Stokes waves  . 3) Breaking of deep-water waves  .
The early nonlinear theories in 2D have now grown into full 3D theories. However, a long-standing problem is that of dimensional reduction of the nonlinear 3D problem to a 2D analysis, as it can be done straightforwardly for the case of shallow water. For deep water waves, the Hamiltonian formulation   provides a nonlinear system for two two-dimensional field variables, namely surface elevation and potential at the surface. To compute the velocity normal to the surface, the solution of the underlying 3D-Laplace equation is achieved by an expansion with respect to powers of the surface elevation into spectral function. These algorithms have been named ‘higher order spectral methods’ (HOSM) and go back to the work of West et al.  and Dommermuth and Yue  . Another attempt is the systematic derivation of an equation describing the slowly varying envelope of a given monochromatic wave train. This equation has the form of a nonlinear cubic Schrödinger equation   and reflects the side-band instability of Stokes waves as well as envelope solitons as localized stable solutions.
Contrary to the Shallow-Water or Thin-Film equations, no intrinsic length scale exists for the deep-water case. The only available scales are given by the solutions sought for, namely the wave length (if the waves are nearly coherent) and the amplitude A of those waves. In deep-water waves, the ratio of amplitude to wave length is normally rather small. Thus, using the wave steepness
with for expanding the nonlinearities seems to be quite natural.
The present paper offers a consistent scheme of reducing the spatial dimension of nonlinear deep-water wave theory by introducing the concept of fractional derivatives along the free surface. Starting from the 2D Euler equations for an incompressible potential flow, a one-dimensional model describing deep-water surface waves is derived. Similar to the Shallow-Water case, the z-dependence of the dependent variables is determined from the Laplace equation and a set of two nonlinear equations for the surface velocity and the surface elevation remains. The system is nonlocal and can be formulated in conservative form, describing waves over an infinitely deep layer. It contains nonlinearities of the lowest (quadratic) order in velocity and surface elevation but can be systematically extended to higher orders. Numerical solutions are presented for different initial conditions. Coherent wave trains prove to be unstable due to the Benjamin-Feir instability  and localized solutions in form of envelope solitons are obtained. Finally, a spatially localized and temporally oscillating surface pressure variation is used to act as a wave-maker. The propagation of different wave ensembles is studied and linear and nonlinear results are compared to show the importance of the nonlinear terms.
The conservation of the total energy can be considered as a quality test for both the derived system and the numerical method applied for its solutions. We show that the energy is conserved within a few percents which give us confidence in our derivation and results.
The model can be extended to the case of a 3D potential flow, resulting in a set of three two-dimensional evolution equations.
2. The Model
2.1. Potential Flow
For the sake of simplicity we restrict the treatment first to two dimensions, . We are looking for solutions of the non-dimensional Euler eqs. for an incompressible fluid in the half space , see Figure 1:
where time is scaled by
and lengths with an arbitrary which will be identified later (sect. 3.2) with the characteristic wave length of the initial condition. We assume that the flow field can be derived from a potential
which then, due to (5) must be a solution of the Laplace eq.
Figure 1. Sketch of the system (2D). A potential flow is considered in the half space .
which also implies
be the horizontal velocity component at the free surface located at .
From (3) we derive
Equation (11) serves as a boundary condition (b.c.) for (9) at . The surface function h is determined by the kinematic b.c.
For the following treatment, it is of advantage to formulate (12) also in conservative form
with the flux
For an infinitely deep layer, the asymptotic b.c.
must hold. Hence the general solution of (9) reads
with c.c. as the complex conjugate.
2.2. Nonlocal Expansion
Considering (11) and (13) as evolution Equations, the original 3D (here 2D) problem is reduced to a 2D (here 1D) system in the horizontal coordinates x,y (here x) only. However, to evaluate the flux (14) one needs to know . Taking (16), Equation (14) reads
To determine the amplitudes from , we evaluate (16) at and expand the exponential function up to a given order :
Next we introduce the fractional differential operator of order :
with the nonlocal operator defined as
for every function via its Fourier transform . We note that (20) is one of the definitions of the fractal Laplacian, here in 1D  . Now Equation (18) can be written as
(in the 3D case, must be replaced by the horizontal Laplacian ).
Using the same technique, S from (17) can be expressed as
leading finally to
where is the operator inverse to and is computed from
Thus, the equations (11) and (13) are closed by the relation (23). Note that is a nonlocal operator only for odd n. For n even it can be expressed by spatial derivatives of order n in position space. Moreover, assuming periodic lateral boundary conditions, it is self-adjoint for all n.
2.3. Second Order Model
In the following we shall restrict ourselves to , resulting in a model of second order in the dependent variables . With
and the self-adjoint operator
Equation (23) turns into
where only terms up to the second order in are included.
In (11), the pressure gradient evaluated at turns into
where is the pressure along the free surface and assumed to be constant (external pressure). To find we evaluate (4) at the surface and include only linear terms in , leading to
with the vertical velocity component
along the surface. is related to by (5) according to
where again only linear terms are considered. Finally, we have
where the first term on the right-hand side corresponds to the hydrostatic pressure, the second one to the dynamic part. Using (33), Equation (11) can be cast into the form
After inversion of the operator on the left-hand side of (34), the complete system up to the second order reads
with from (28).
2.3.2. Linear Waves
Taking only linear terms into account one finds from ((35), (36))
or, eliminating , the linear wave equation
where we have used the identity . Equation (38) possesses the well-known short-wave dispersion relation
2.3.3. Two Horizontal Dimensions
If we add the second horizontal dimension we are able to describe a 3D flow which, however, is still potential. Extending the surface velocity to a 2D vector
the derivation of the reduced system is straightforward and it reads
with the 2D nabla operator .
2.4. Similarities and Differences with HOSM
Contrary to the higher-order spectral methods (HOSM) presented in   , we use the surface elevation and the surface velocity instead of the surface potential. With the kinematic b.c. (36) formulated by the help of the flux S, this results in model equations that can be written in conservative form. Due to the self-adjointness of it can be easily seen that are conserved in time, also for the numerical scheme, if the derivatives with respect to x are expressed by simple central differences (see sect. 3.1 below). On the other hand, our method can be formulated in a closed form in position space and is, at least up to second order, faster and easier to implement numerically, allowing for larger aspect ratios and an effective numerical realization in two horizontal dimensions. Eventually, Equations ((35), (36)) can be considered as a systematic expansion with respect to the wave steepness . With the self-consistent assumption , our model includes all terms up to the order . Hence it should deliver accurate results up to second-order stokes waves, which corresponds to  , about 36 percents of the steepest possible Stokes wave with . Then the maximal wave amplitude of a wave with (wave length ) should not be much higher than .
3. Numerical Method and Results
3.1. Numerical Method
We concentrate on the 2nd order model ((35), (36)) with x in the domain and periodic lateral b. c. . To evaluate (28), the expressions
must be computed. The operator introduced in (20) thus reads . It is nonlocal for odd n and (45) can be written as
with the Green’s function
and . Taking a finite difference scheme for the numerical integration of the system, the variable x is discretized as
where N is the number of mesh points and with . The discrete form of (46) thus reads
which would result in time-consuming computations of convolutions when evaluating (28). Hence it is much more effective to compute (48) in Fourier space where it simply reads
where is the discrete Fourier transform of that can be obtained numerically by any standard Fast-Fourier transform (FFT)  .
To evaluate S according to (28) it is therefore necessary to compute and in Fourier space, then perform the products and in real space and finally compute in Fourier space again. In addition for (35), has to be determined in Fourier space. Thus, seven FFTs (three forward and four backward) have to be performed at each time step.
The derivatives with respect to x occurring in ((35), (36)) are approximated with centered differences
what ensures the conservation of .
Finally, time integration is achieved by an explicit 4th-order Runge-Kutta algorithm  .
3.2.1. Benjamin-Feir Instability of Slightly Disturbed Coherent Wave Trains
The initial conditions are based on surface elevations in form of a traveling wave with given amplitude and wave number . If for the spatial scaling introduced in (6) the wave length is chosen, we can put . The x-dimension is discretized with mesh points, the step sizes used are , , resulting in a side length of and allowing for 60 initial wave lengths. Each wave length is resolved with about 33 mesh points.
As first shown in  , coherent Stokes waves with a single frequency and wave length are unstable with respect to the side-band instability. The time series in Figure 2 presents the evolution of a 2nd-order Stokes wave  traveling to the right side as initial condition:
Figure 2. Snapshots of a BF unstable wave train. Numerical solution of the 1D model ((35), (36)).
To trigger the side band (Benjamin-Feir, BF) instability  , weak disturbances of the form
with and have been added. The speed and wave length of the (undisturbed) traveling wave follows as
what yields a Courant number of
The disturbances (39) belong to neighbored waves of in Fourier space. The absolute values of the Fourier amplitudes over time for the carrier wave ( ) as well as for its five neighbors ( ) are shown in Figure 3. In the beginning, the side bands grow exponentially, taking the energy from the carrier wave that decreases.
3.2.2. Envelope Solitons
If a coherent wave is enveloped by a slowly varying function according to
Zakharov  showed that for small but finite amplitudes an equation for the complex valued can be derived, having the form of the nonlinear Schrödinger equation (NLS):
Figure 3. Absolute values of the Fourier amplitudes over time for the run shown in Figure 2. Solid: carrier wave , dashed-dotted: side bands , dashed: side bands .
The nonlinear coefficient depends on the ratio of the water depth to the wave length of the carrier wave  and reads for infinite depth
Equation (56) has the well-known localized hyperbolic secant solutions
with the relation between reciprocal width and amplitude
To see how envelope solitons behave in our model, we take as initial condition
Figure 4 shows snapshots from the evolution for
with and . The envelope soliton travels thereby with the group speed as given in (57) which is half of the phase speed of the carrier wave.
To confirm the agreement of our model and its numerical solutions with the results of the weakly nonlinear NLS, we also studied an initial condition as (61), but with an amplitude and width of the envelope that does not satisfy the relation (60). Taking the same but , the package flows apart and finally extends over the whole domain, Figure 5.
Finally we show in Figure 6 the interaction of two colliding envelope solitons, generated by the initial condition
The solitons travel in opposite direction and interact without changing their shapes strongly. Due to the periodic b.c. these interactions repeat after equidistant time intervals of .
3.2.3. Conservation of Energy
It is interesting to check the temporal variation of the energy during the runs. For an inviscid fluid, the total energy
Figure 4. Envelope soliton for an initial condition (61) fulfilling the relation (60).
with the potential energy
and the kinetic energy
must be conserved. For the latter one, one finds up to quadratic order the expression
Figure 5. Evolution of a localized initial condition (61) which does not fulfill the relation (60) for envelope solitons.
Figure 6. Two envelope solitons with different amplitudes traveling in opposite directions.
The conservation of (63) can serve as a quality feature for the numerical solutions as well as for the accuracy of the derived system itself.
Figure 7, Figure 8 show the three energies (63)-(65) over time for the series depicted in Figure 2 and Figure 6, respectively. In Figure 8, potential as well as kinetic energy varies strongly during the collision phases but the sum (63) remains fairly constant. Table 1 shows the standard deviations normalized with the square mean defined as
Figure 7. Energies plotted over time of the evolution shown in Figure 2. solid: total, dashed: kinetic, dashed-dotted: potential.
Figure 8. Energies of the evolution shown in Figure 6. solid: total, dashed: kinetic, dashed-dotted: potential.
Table 1. Standard deviations of the energies shown in Figure 7, Figure 8.
with the mean
and f as one of the three energy functions.
Note that all computations are only up to second order which then applies also for the energy conservation so that
where stands for h, and combinations.
3.2.4. High Waves Caused by Dispersive Focusing of an Unstable Stokes Wave
Waves on deep water are generated by wind and have a broad spectrum of wave lengths. Since their phase and group speed is proportional to , longer waves are faster than shorter ones and, if they propagate in the same direction, may catch up and enhance the total amplitude. This is already a linear effect but can be amplified by nonlinearities what then leads to even higher wave heights . In this way, very large heights can be obtained, which is one of the mechanisms responsible for the occurrence of freak or rogue waves  and named “dispersive focusing” or “dispersive enhancement”.
To classify waves with respect to their height, the “abnormality index” AI
is introduced, where is the height of the wave and the significant wave height. The significant wave height is defined as the average height of the upper third of the wave height distribution (over time or space). Another way to compute is to use the root of the zeroth order moment of the spectrum
Then, a freak wave can be defined by AI > 2. The probability for freak waves depends strongly on the initial wave pattern and on the dimensionality of the system. In 2D, only two propagation directions are possible and freak waves should occur quite often. Figure 9 shows the evolution of an unstable Stokes wave as (52), but for a larger aspect ratio and a higher grid resolution at . A freak wave shows up with AI ≈ 2.8. In Figure 10 we plot AI over time, the circle marks the situation shown in Figure 9.
with can be seen over wide regions, in agreement with the findings of Zakharov et al.  .
3.2.5. Waves Generated by Localized Pressure Forcing
In this section we consider an initially quiescent flat layer with at
Figure 9. A large wave with AI ≈ 2.8 is formed after from the initial condition (52). Here, and , resulting in a spatial resolution of almost 70 mesh points per carrier wave length.
Figure 10. Abnormality index over time for the run leading to Figure 9, computed along (70) with (71). Two singular events are evident, the other one around .
Figure 11. Absolute values of the Fourier amplitudes of shown in Figure 9.
on which waves are generated by a spatially localized surface pressure of the form
Inclusion of (73) modifies (35) to
If the amplitude varies harmonically, (73) acts as a wave maker and produces waves with wave length , traveling symmetrically to both sides. To simulate the situation where waves with shorter wave lengths are followed by longer (faster) waves, we perform a run with and . Frequency and amplitude of the pressure are decreasing linearly in time, according to
with . Two runs have been performed, one with the nonlinear model ((74), (36)), the other with the same system but only keeping linear terms in . From Figure 12 it is evident that the nonlinearities amplify the amplitudes and provide significantly higher waves. Maximal focusing is obtained after at , see Figure 13.
We presented the derivation and numerical solutions of dimension-reduced deep-water equations. The system can be considered as a systematic expansion
Figure 12. Abnormality index of the wave field over time for waves generated by a local pressure variation (73) with frequency and amplitude variation according to (75). Black: nonlinear model (2nd order), blue: linear equations. Amplitudes are significantly enhanced by nonlinear coupling. Maximal focusing takes place after at .
Figure 13. Snapshots before, while and after maximal focusing. mesh points are used for , but only half of the layer is shown.
of the basic hydrodynamic equations with respect to powers of the wave steepness. We discussed the model in the lowest nonlinear order and showed numerical 1D-solutions where the Benjamin-Feir instability of coherent waves turned out and the behavior of envelope solitons has been studied. All results are in good agreement with previous work. The envelope solitons coincide well with solutions of the nonlinear Schrödinger equation both in amplitude and width. Our model is not based on a certain carrier wave and no spatial direction along the surface is distinguished. The formulation is rotationally invariant in the horizontal dimension(s) and allows therefore also for the study of colliding wave trains.
From the numerical side, the great advantage of the present model is its relatively straightforward and effective numerical realization. Due to the dimension reduction, CPU-times of the runs shown are in the range of hours on normal PCs or Laptops. Our system can be extended both to higher order nonlinearities as well as to three-dimensional potential flows, resulting in a reduced 2D set. Thus, the numerical examination of 2D envelope solitons occurring on the surface of a 3D potential flow should become feasible. This work is currently in progress and will be presented in a following paper.
 Saint-Venant, A.J.C. (1871) Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l’introduction de marées dans leurs lits. Comptes Rendus des Séances de Académie des Sciences, 73, 147, 237.
 Bestehorn, M. (2009) Fluid Dynamics and Pattern Formation. In: Meyers, R.A., Ed., Contribution to Encyclopedia of Complexity and System Science, 4, Springer, Berlin, 3611.
 Longuet-Higgins, M.S. and Cokelet, E.D. (1976) The Deformation of Steep Surface Waves on Water. I. A Numerical Method of Computation. Proceedings of the Royal Society of London, A350, 1.
 Zakharov, V.E. (1968) Stability of Periodic Waves of Finite Amplitude on the Surface of a Deep Fluid. Journal of Applied Mechanics and Technical Physics, 9, 190.
 West, B.J., Brueckner, K.A., Janda, R.S., Milder, D.M. and Milton, R.L. (1987) A New Numerical Method for Surface Hydrodynamics. Journal of Geophysical Research, 92, 11803.
 Zakharov, V.E., Dyachenko, A.I. and Prokofiev, A.O. (2006) Freak Waves as Nonlinear Stage of the Stokes Wave Modulation Instability. European Journal of Mechanics, 25, 677.