Theory of deep-water surface gravity waves derived from a Lagrangian

An exact set of equations describing deep-water irrotational surface gravity waves, originally proposed by Balk (Phys. Fluids, vol. 8 (2), 1996, pp. 416–420), and studied in the case of standing waves by Longuet-Higgins (J. Fluid Mech., vol. 423, 2000, pp. 275–291) and Longuet-Higgins (Proc. R. Soc. Lond. A, vol. 457 (2006), 2001, pp. 495–510), are analytically examined and put in a form more suitable for practical applications. The utility of this approach is its simplicity. The Lagrangian is a low-order polynomial in the Fourier coefficients, leading to equations of motion that are correspondingly of low degree. The structure of these equations is examined, and the existence of solutions is considered. To gain intuition about the system of equations, a truncated model is first examined. Linear stability analysis is performed, and the evolution of the fully nonlinear system is discussed. The theory is then applied to fully resolved permanent progressive deep-water waves and a simple method for finding the eigenvalues and eigenvectors of the normal modes of this system is presented. Potential applications of these results are then discussed.


Introduction
Nonlinear surface wave models are a powerful tool for studying complex wave scenarios, and in particular for investigating the properties of steep and overturning waves. An improved description of the dynamics of surface water waves, and their behaviour under these nonlinear conditions, is crucial for a better understanding of airsea interaction processes (Melville 1996). Nonlinear surface waves pose a formidable theoretical problem, and analytical studies of weakly nonlinear waves have led to a better understanding of water wave phenomena, but these results are limited by their restriction to weak nonlinearities. That is, many of these theoretical predictions are far outside of their region of validity as wave breaking is approached and during the subsequent breaking event. For this reason, numerical models of surface waves can provide insight into the physics of highly nonlinear waves processes. In this paper we derive, and put in a form more suitable for numerical applications, the model of Balk (1996), based on a Lagrangian for deep-water surface gravity waves.
The outline of this paper is as follows. First, in § 2 we provide a derivation of the model of Balk. In § 3 the equations are examined. Next, in § 4 we look at a truncated model, with only one spatial mode. In § 5, the linear stability of Stokes waves is examined. The results are then discussed in § 6.

Derivation of the equations of Balk
In this section we derive the governing equations for irrotational inviscid surface gravity waves in a conformally mapped reference frame. This formulation has a number of distinct advantages, including the ability to model wave overturning, a knowledge of the truncated properties (a necessary feature of any numerical model) of the waves (e.g. energy, momentum), and a relatively straightforward recipe for numerical implementation. Although the model is algebraically intensive, the tools needed to derive these results are relatively simple (cf. the boundary integral methods of Longuet-Higgins & Cokelet (1976) and Dold & Peregrine (1986)). Instead of following Balk directly, we make use of some results from Dyachenko et al. (1996) and show how this approach leads to a more natural derivation of the governing equations originally found by Balk. To this end, we consider irrotational two-dimensional inviscid flow in a fluid of infinite depth with a free surface. Density and the gravity constant are set to 1 and the flow is assumed to be 2π periodic in the x-coordinate.
The domain in the complex plane of the fluid is given by z = x + iy, which can be considered to be the conformal image of the lower half-plane described by ζ = ξ + iκ (these variables are time dependent, but this is not explicitly written now for clarity of presentation).
The transformation from the (ξ , κ) plane to the (x, y) plane is given by z(ζ ) = x(ξ , κ) + iy(ξ , κ). The velocity of fluid w = u + iv is given by the complex potential f (z) = ϕ(x, y) + iΨ (x, y) (here ϕ and Ψ represent the velocity potential and streamfunction, respectively, and we assume they satisfy the condition v → 0 as y → −∞). That is, where C(ζ ) = a(ξ , κ) + ib(ξ , κ) is the complex potential in the ζ plane and * denotes the complex conjugate. The surface of the fluid is going to be the image of the real axis in the ζ -plane, i.e. 2a,b) and the surface values of the velocity potential and streamfunction are We note that {X, Y, A, B} completely determine the state of the fluid. Furthermore, Balk notes that knowledge of {Y(ξ ), B(ξ )} are sufficient to completely characterize the fluid, since the boundary values at κ = 0 of the real parts of analytic functions in the lower half-plane can be found from the imaginary parts via the Hilbert transform.
Recall, in the x-y plane, we have from Zakharov (1968, see also Miles 1977) that the action for water waves can be written as

N. Pizzo
Now, from the chain rule the time derivative of the free surface η can be rewritten in terms of the conformal variables as (cf. Dyachenko et al. 1996, equation (2.14)) where dots and primes over the summation symbol represent differentiation with respect to time and ξ , respectively. Therefore, this term in the action becomes Next, the kinetic and potential energy densities, per unit wavelength, can be written as while the kinetic energy density is given by (via a Green integral identity) Finally, exploiting the analyticity of C, we can use the Cauchy-Riemann equations to rewrite (2.4) as with the time interval (t a , t b ) chosen so that the dependent variables are fixed at these times. Now, variations of S with respect to A give which is precisely the kinematic boundary condition, in the conformally mapped coordinates (equation (2) in the paper of Balk). The constraint that mass is conserved is the condition so that taking the constant of integration to be 0, corresponding to the mean water level being at y = 0, we find 1 2π 2π 0 YX dξ = 0, (2.12) which is equation (3) in Balk. Substituting (2.10) into (2.9) we find that the action becomes with the constraint given in (2.12). In the derivation given above, the constraints on the system enter naturally by a change of variables in the action, whereas in Balk's derivation, he simply imposes these conditions. Note, a priori there is no reason that the Lagrangian should be the difference of kinetic and potential energies (e.g. see (2.4)) in an Eulerian framework, which has led to significant confusion in the literature (see, for example, the discussion in Luke (1967) and Salmon (1988)). However, both Balk's original formulation and equation (8) in Luke (1967) constrain the system to obey conservation of mass, the kinematic boundary condition, and the condition of no flow at the bottom. When these conditions are met, then the Lagrangian for water waves may be set equal to the difference of kinetic and potential energies (see the discussion in Luke (1967)).
Conservation of mass can explicitly be built into the Lagrangian by exploiting the periodicity of the flow in the x-direction. This property implies the following Fourier expansions: (2.14) and where σ k = (1, 0, −1) for (k > 0, k = 0, k < 0), respectively. This gives us relationships between the real and imaginary parts of the conjugate pairs. We note explicitly where Re/Im stand for the real and imaginary parts, respectively. Using the kinematic boundary condition, i.e. equation (2.10), we find (for n = 0) The complex potential is defined up to an arbitrary constant, so that without loss of generality, we may set A 0 = B 0 = 0. Substituting the Fourier expansions of {A, B} into the equation for the kinetic energy, i.e. equation (2.8), we find Similarly, with these expansions the potential energy is found directly from (2.7) and (2.14) (see also Longuet-Higgins 2000) The Lagrangian, L = T − V, can now be written as where B k is given in (2.20) and we choose Y 0 to satisfy the constraint given in (2.12), namely, The two most obvious conserved quantities associated with this Lagrangian are the x-momentum (associated with phase shift invariance), defined as and energy (associated with time shift invariance) given by Next, the governing equations are found by varying the Lagrangian (2.23) with respect to Y k . This yields the Euler-Lagrange equations for each mode k d dt For partial wave solutions, the infinite sums used in the Fourier expansion of the dependent variables will be truncated at a finite N, so that we will have a truncated Lagrangian L N = T N − V N . This is advantageous since we know the conserved quantities associated with this truncated Lagrangian before deriving the dynamical evolution equations (see Salmon (1988) for a thorough discussion).

Equations of motion
In order to better understand the governing equations, we extend the work of Longuet-Higgins (2000, 2001 from the special case of standing waves to the more general equations of motion. First, from (2.21) we have for n > 0, and similarly for n < 0, (3.2) By substituting these two expressions into (2.21) we find that we can write the kinetic energy as The operator [ * ] is defined as and Y k ≡ 0 when k > N for N the number of Fourier modes in our model, i.e. the resolution. Note, by construction, the first term in parentheses in (3.3) is the complex conjugate of the second term. By using these definitions, and reversing the order of summation, the kinetic energy can be rewritten as That is, the matrix Q kl is related to the column by column multiplication of two P matrices. We now use these definitions to solve for the equations of motion. Substituting the kinetic energy into the first term of equation (2.27)  (3.10) Putting this together, we see that our governing equations are now of the form (3.11) for each n = (±1, ±2, . . . ± N).
Further simplifications of these equations can be made by realizing that ( 3.12) From the definition of P ij , i.e. equation (3.4), we find that Next, we consider the potential energy term, which, from equation (2.7), can be written as where F(i, j, k) = |ijk| : i + j + k = 0 and i, j, k = 0, 0 : otherwise.
Finally, we define where and Einstein summation is assumed. Together with initial conditions at time t = 0 for Y n andẎ n , equation (3.16) are the complete equations for surface gravity waves written in a compact form. Analytically, we may also check that our equations reduce to those of Longuet-Higgins (2000), for the case where all Y k are real, which corresponds to standing waves. It is a simple matter to show that under the condition that Y −k = Y k , and by mapping into the coordinates used in that work, namely a n = 2|n|Y n with a 0 = 1, then the entries of P become P mn → 1 2 1 √ |m||n| (a n−m + a m+n − a m a n ).
( 3.18) With this identification, the rest of our variables can easily be shown to be equivalent to those of Longuet-Higgins (2000).

Example: N = 1
In this section we consider a simple example where only one mode is present. We perform linear stability analysis of permanent progressive waves, as this system illustrates many properties of the large N Stokes waves to be considered in the next section. Furthermore, we examine properties of general solutions to this system of equations.

N. Pizzo
Note, these governing equations can be rewritten as where F 1 , F 2 are the right-hand sides of equations (4.3) and (4.4), respectively. We define the coefficient matrix Q as (4.10) Now, equation (4.9) has solutions provided that Q is invertible. Therefore, there is some interest in looking at how the wave evolution breaks down as det(Q) → 0. To this end, the determinant of Q is The reality condition on the Fourier components mandates that α * = β, so that letting α = A + iB = β * , this condition implies that (A, B) must lie within a circle of radius 1/2, i.e.

Permanent progressive waves
A simplified class of solutions to equation (4.9) are those of the permanent progressive type. That is, we choose (4.13a,b) where c is a to be determined phase velocity of the waves. When this ansatz is made, we find from (4.9) the following relationship between c and α 0 : (4.14) For α 0 small we approach the infinitesimal sinusoidal wave solution we expect from linear theory for deep-water surface gravity waves (Phillips 1977). As α 0 → 1/2, the free surface forms a cusp (note this is also when the coefficient matrix Q becomes singular), taking the form of a cycloid (see figure 1). This may also be seen as the point where X ξ (ξ = π) = Y ξ (ξ = π) = 0.
The total energy of these waves T + V is conserved, and takes the form Similar to Stokes waves (Schwartz 1974;Longuet-Higgins 1975), the energy does not monotonically increase with α 0 . The energy maximum occurs when α 0 = 1/ √ 6. This has implications for the stability of the system, as is discussed below.

Normal mode stability analysis
The linear stability analysis of this system is performed by letting where 1 is an ordering parameter. Applying Hamilton's principle to the Lagrangian at second order in , i.e. O( 2 ), we find where (M, K, C) are 2 × 2 matrices and A = (A, A * ). Defining ν = −4α 2 0 and µ = ν(1 − 2α 2 0 ), the matrices are given by Note, M and C are symmetric while K is skew-symmetric. We now let A = A 0 e λt , where λ is the normal mode eigenvalue and (A 0 , A * 0 ) the corresponding eigenvector. The eigenvalues are found by requiring that the determinant of the left-hand side of equation (4.17) vanishes. That is This implies that the eigenvalues are given by (4.20) We see immediately that the eigenvalues are real when α 0 > 1/6, (4.21) which is equivalent to the point at which the energy becomes a maximum. This is in agreement with existing studies on Stokes waves (Saffman 1985;MacKay & Saffman 1986). The corresponding eigenvectors are found by solving equation (4.17) with the eigenvalues substituted in. These are unique up to a rescaling by a complex constant. When Λ < 0, solutions may be unstable, leading to the formation of a cusp in the free surface as ρ → 1/2.
For ρ small we see the potential goes to positive infinity, as P 2 0. For ρ → 1/2, the sign of the potential is given by the sign of Λ = P 2 − H + 1/8. (4.28) When Λ < 0, solutions might lead to cusp formation as the value of ρ can approach 1/2.

Stokes waves
The normal mode linear stability of Stokes waves has been studied in detail by, for example, Longuet-Higgins (1978a), McLean et al. (1981), Tanaka (1983) and MacKay & Saffman (1986). Here, we present an alternative derivation of the linear stability analysis with the added benefit that the system of equations takes a particularly simple form and is general. That is, once the Stokes coefficients (this may, for instance, be for class M Stokes waves (Saffman 1980)) are determined, the stability analysis follows readily by solving an eigenvalue problem with prescribed matrix entries which are low-order polynomials in the dependent variables, making this approach particularly suitable for numerical implementation.

Stokes waves
We now consider permanent progressive solutions to the system where N is taken to large enough to properly resolve the physical scenario in question.
To this end, we note that Stokes wave solutions to our system of equations take the simple form (Balk 1996) Y k = α k e ickt , (5.1) so that B k = cY k . The Lagrangian can then be written as with V given by equation (2.7). Therefore, the system governing the coefficients α k is (Longuet-Higgins (1985), see also Longuet-Higgins (1978c)) for F n = 0. There are N + 1 equations here, so we must prescribe a parameter to close the system. Generally, the phase speed c has been used. However, c does not increase monotonically with increasing slope and contains critical points, which leads to loss in numerical accuracy (see § 5 of Longuet-Higgins (1985)). Alternatively, a variety of other parameters which do increase monotonically with increasing ak have been proposed. Following Longuet-Higgins (1985), we take the parameter Q = 1 − q 2 crest , where q crest is the particle speed at the crest. From Bernoulli's equation, we have Q = 1 + Y(0) = 1 + 1 2 α 0 + α 1 + α 2 + · · · . To solve the system of equations, i.e. (5.3), we employ a Newton-Raphson-type method.

Linear stability analysis
To find the form of the linear normal modes, we next let where is a small parameter, λ is the normal mode eigenvalue and A k are the coefficients of the eigenvector.
The Lagrangian to O( 2 ) becomes and primes represent sums over all indices except 0. Explicitly, the mean surface displacement is given by so that its time derivative iṡ Substituting these relationships into the B k term, we find where superscripts imply we are taking terms of that order from the expansion. This is sufficient to define the Lagrangian. Applying the Euler-Lagrange equations, we find that the governing equation takes the form (cf. the discussion for the N = 1 case in § 4) of a quadratic eigenvalue problem (Tisseur & Meerbergen 2001). Namely, we have (λ 2 M + λC + K)A = 0,  where G(n, m) = 1 if −N − m n N − m and 0 otherwise. The entries of these matrices are low-order polynomials in the Stokes coefficients, making them particularly suitable for numerical implementation. Furthermore, these results are general in the sense that one may examine subharmonic instabilities (Longuet-Higgins 1978b) by considering class M, for M a positive integer, Stokes waves (Saffman 1980;Zufiria 1987).
To validate these algebraically complex computations, the eigenvalues corresponding to superharmonic perturbations of steep Stokes waves are computed from this method and are shown on top of the previous computations due to Tanaka (1983) in figure 3. Agreement between the predictions of equation (5.12) and Tanaka's results are found. Furthermore, the change in (minus) total energy versus slope, ak, is shown in the plot, showing that the change in stability corresponds to the critical point in the energy.

Conclusion
In this paper we have transformed the equations of Balk into a form more suitable for theoretical and numerical examination. Furthermore, we have derived the equations of motion that result from Balk's Lagrangian. This represents a generalization of the work done by Longuet-Higgins (2000, 2001 for standing waves. The system of equations are found for independent coefficients Y n (t), n = 1, 2 . . . , with partial wave solutions found by truncating the system at n = N. Solutions exist as long as the matrix Q, i.e. equation (3.17), is invertible.
A one mode example was considered, and the stability of permanent progressive waves was examined. Linear stability conditions were derived, and the evolution of the fully nonlinear system was considered. The linear stability of Stokes waves (here the number of modes N was increased until convergence was found) was then inspected. A general system of equations to derive the eigenvalues was presented, and comparison was made to the values found by Tanaka (1983).