\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,published]{cedram}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\usepackage{etoolbox}
\usepackage{multirow}
\usepackage{upgreek}
\csdef{Seqnsplit}{\\}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20230171}
\def\morerows#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\def\botline{\\\hline}
\def\balpha{\boldsymbol{\alpha}}
\def\bfbeta{\boldsymbol{\beta}}
\def\bvarepsilon{\boldsymbol{\varepsilon}}
\def\bsigma{\boldsymbol{\sigma}}
\def\unient#1{}
\def\HCode#1{}
\def\mathbi#1{\textbf{\textit{#1}}}
\allowdisplaybreaks
\DOI{10.5802/crmeca.221}
\datereceived{2023-02-15}
\daterevised{2023-07-17}
\dateaccepted{2023-07-13}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage{multirow}
\usepackage{stmaryrd}
\def\CDRdoi{10.5802/crmeca.164}
\def\0{\phantom{0}}
\def\newline{\unskip\break}
\def\sfrac#1#2{{#1}/{#2}}
}

\makeatother

\dateposted{2023-08-03}
\begin{document}


\begin{noXML}

\TopicEF{Research article}{Article de recherche}

\title{Interfacial capillary--gravity short-crested waves}

\alttitle{Les ondes \`a courtes cr\^etes interfaciales de
gravit\'e--capillarit\'e}

\author{\firstname{Dalila} \lastname{Boughazi}}
\address{Universit\'{e} M'Hamed Bougara, Facult\'{e} des Sciences,
Si\`{e}ge (Ex. INIL) Boumerd\`{e}s 35000, Alg\'{e}rie}
\address{USTHB, Facult\'{e} de physique, BP 32 El Alia Alger, Alg\'{e}rie}
\email[D. Boughazi]{bouda.2007@hotmail.fr}

\author{\firstname{Mohammed} \lastname{Debiane}\IsCorresp} 
\addressSameAs{2}{USTHB, Facult\'{e} de physique, BP 32 El Alia Alger, Alg\'{e}rie}
\email[]{debianemd@yahoo.fr}
\email[M. Debiane]{mddebiane@gmail.com}

\author{\firstname{Nabil} \lastname{Allalou}}
\addressSameAs{2}{USTHB, Facult\'{e} de physique, BP 32 El Alia Alger
Alg\'{e}rie}
\email[N. Allalou]{n\_allalou2004@yahoo.fr}

\keywords{\kwd{Short-crested interfacial waves}
\kwd{Gravity}
\kwd{Capillarity}
\kwd{Arbitrary depths}
\kwd{Resonance}}

\altkeywords{\kwd{Ondes \`a courtes cr\^etes interfaciales}
\kwd{Gravit\'e}
\kwd{Capillarit\'e}
\kwd{\'Epaisseurs arbitraires}
\kwd{R\'esonance}}

\begin{abstract}
This paper presents a method derived from Whitham's variational
formulation of the problem of interfacial capillary--gravity
short-crested waves. It is developed for the resolution of the problem
of waves generated by obliquely reflecting interfacial waves from a
vertical wall. In essence, Whitham's method is not changed, but
computations are performed and arranged to produce a method that has
been applied to a number of cases demonstrating the effectiveness and
the flexibility of the approach. The performance of the method is
illustrated in several examples including
the case of harmonic resonance.
\end{abstract}

\begin{altabstract}
Cet article pr\'esente une m\'ethode d\'eriv\'ee de la formulation
variationnelle de Whitham du probl\`eme des ondes interfaciales \`a courte
cr\^ete de gravit\'e--capillarit\'e. Elle est d\'evelopp\'ee pour la
r\'esolution du probl\`eme des ondes g\'en\'er\'ees par la r\'eflexion
oblique des ondes interfaciales sur une paroi verticale. La m\'ethode
de Whitham n'est pas modifi\'ee, mais des calculs sont effectu\'es et
organis\'es pour produire une m\'ethode qui a \'et\'e appliqu\'ee \`a un
certain nombre de cas d\'emontrant l'efficacit\'e et la flexibilit\'e
de l'approche. La performance de la m\'ethode est illustr\'ee par
plusieurs exemples, y compris le cas de la r\'esonance 
harmonique.\looseness=-1
\end{altabstract}

\shortrunauthors

\maketitle
\vspace*{-.8pc}
\end{noXML}

\section{Introduction}\label{sec1}
Over the past decades free surface short crested waves have been
studied extensively. These waves are one of the simplest classes of
three-dimensional water waves and which are periodic 
 in each of two
distinct horizontal directions. They arise at the linear level from the
superposition of two progressive wave-trains of equal amplitude and
frequency which are propagating at an angle to each other. Starting
with the work of Fuchs~\cite{1} who obtained second-order solutions,
approximations through perturbation expansions of higher orders were
performed, i.e.,  by  Hsu \etal~\cite{2},  Roberts~\cite{3}
Ioualalen~\cite{4}, Kimmoun \etal~\cite{5},  Craig and
Nicholls~\cite{6,7}. Roberts and Schwartz~\cite{8} and
Okamura~\cite{9,10} computed fully numerical short-crested wave
solutions. Debiane and Kharif~\cite{11} calculated resonant and
non-resonant short-crested gravity--capillary waves in deep water using
a method derived from Whitham's variational formulation.

While free surface short-crested waves have been extensively studied,
relatively little theoretical works has been made for internal
short-crested waves. To our knowledge, only two studies have been
devoted to this field. Allalou \etal~\cite{12} have investigated
numerically fully  three-dimensional gravity waves at the interface
between two fluids of arbitrary depths confined between two horizontal
rigid walls. They computed solutions via a perturbation expansion up to
the 27th order in wave steepness. Bocharov \etal~\cite{13} have
numerically solved the weakly nonlinear model equation proposed by
Khabakhpashev and Tsvelodub~\cite{14}. The shapes of stationary
traveling internal waves, both periodic in the two horizontal
coordinates and soliton-like, are presented.

Motivated by the efficiency of the method proposed by Debiane and
Kharif~\cite{11} to the problem of surface short crested waves, we
develop here a similar method for the problem of short-crested waves
on the interface of a two-layer fluid, in the presence of gravity and
surface tension. Even though density stratifications that occur in
nature are continuous, it is reasonable to use the two-layer fluid
model with a density discontinuity due to its simplicity and the fact
that it provides a good approximation when the interfacial wavelength
is sufficiently longer than the length scale of the density variation.
The study of these waves is the first step to develop an understanding
of understanding of little investigated three-dimensional interfacial
water waves.

\section{The governing equations}\label{sec2}

We consider short-crested wave field on the interface of a two-layer
system of fluids bounded by two rigid horizontal walls. These waves are
generated by an interfacial wavetrain of wavelength ${\lambda}$
which arrives at some angle of incidence $\theta $ on a vertical
wall and which is perfectly reflected (Figure~\ref{fig1}). $\theta$  is the
angle between the direction of propagation of the incident wave and the
normal to the wall. The resulting waves propagate with constant 
speed $c$ along the wall. The wavelength in the $x$-direction is $\lambda
_{x}=\lambda /\sin \theta $ and the wavelength in the $y$-direction is
$\lambda _{y}=\lambda /\cos \theta $. The fluids are supposed to be
homogeneous, incompressible and inviscid, and the motion is assumed to
be irrotational. The physical variables associated with the lower fluid
are denoted by a subscript~(\ref{eq1}) and those of the upper fluid by a
subscript~(\ref{eq2}). A cartesian coordinate system $\Re
(o,x,y,z)$ is adopted with the $x$- and $y$-axis located on the horizontal
plane, and the $z$-axis pointing vertically upwards. The level 
$z = 0$ coincides with the separation surface between the two layers at rest.
All the variables and the equations will be written in dimensionless
form in which ${K}^{-1}=\lambda /2\uppi$ and
$({g}{K})^{-\sfrac{1}{2}}$ are the length and time of
reference with $K$ being the wavenumber. It is computationally convenient
to use the frame of reference $\Re '(O,X,Y,Z)$, where
{\begin{equation}\label{eq1}
\left(\begin{array}{@{}c@{}}
{X}\\
{Y}\\
{Z}
\end{array}\right)=\left(\begin{array}{@{}c@{}}
{px}-{\varpi t}\\ {qy}\\ {z}
\end{array}\right)
\end{equation}}\unskip
and in which the short-crested wave is steady with a period of $2\uppi$
in $X$ and $Y$. Here $\varpi$ is the non-dimensional frequency and $p$
and $q$ are the non-dimensional wavenumbers in the $X$ and $Y$
directions, respectively, defined by
{\begin{equation*}
p=\sin\theta \quad\text{and}\quad q=\cos \,\theta.
\end{equation*}}\unskip

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig1}Sketch of the reflection of plane wave onto
vertical wall at angle ${\theta}$ of incidence.}
\end{figure}


With the assumptions of incompressibility and irrotational flow in each
of the layers, the fluid motions can be described by velocity
potentials $\Phi _{1}({X},{Y},{Z})$ and $\Phi _{2}({X},{Y},{Z})$ 
in the lower and the upper layers, respectively. These potentials
satisfy Laplace's equation in the lower and upper fluid domains
{\begin{equation}\label{eq2}
{p}^{{2}}\Phi _{{iXX}}+{q}^{{2}}\Phi _{{iYY}}+
\Phi _{{iZZ}}={0}\quad ({i}={1,2}).
\end{equation}}\unskip

These equations are subject to boundary conditions on the interface
between the layers defined by the equation ${Z}=\eta ({X},{Y})$
{\begin{eqnarray}\label{eq3}
&&-\varpi \eta _{{X}}+{p}^{2}
\eta _{{X}}\Phi _{{iX}}+{q}^{2}
\eta _{{Y}}\Phi _{{iY}}-\Phi _{{iZ}}=0 \Seqnsplit
\label{eq4}
&&-\frac{\varpi }{\mu -1}(\mu \Phi _{2{X}}-\Phi _{1{X}})+
\frac{\mu }{2(\mu -1)}({p}^{2}{\Phi }_{2{X}}^{2}+
{q}^{2}{\Phi }_{2{Y}}^{2}+{\Phi }_{2{Z}}^{2})\nonumber\\
&&\qquad -\,\frac{1}{2(\mu -1)}\left({p}^{2}{\Phi }_{1{X}}^{2}+
{q}^{2}{\Phi }_{1{Y}}^{2}+{\Phi }_{1{Z}}^{2}
\right)+\eta -\kappa \xi ({X},{Y})=0
\end{eqnarray}}\unskip
where 
{\begin{eqnarray}\label{eq5}
&\displaystyle\mu =\rho _{2}/\rho _{1} \Seqnsplit
\label{eq6}
&\displaystyle\kappa =\frac{{TK}^{2}}{{g}
(\rho _{1}-\rho _{2})}\Seqnsplit
\label{eq7}
&\displaystyle\xi ({X},{Y})=\frac{{p}^{2}\eta _{{xx}}
(1+{q}^{2}{\eta }_{{Y}}^{2})+{q}^{2}
\eta _{{YY}}(1+{p}^{2}{\eta }_{{X}}^{2})-
2{p}^{2}{q}^{2}\eta _{{XY}}\eta _{{X}}
\eta _{{Y}}}{(1+{p}^{2}{\eta }_{{X}}^{2}+{q}^{2}
{\eta }_{{Y}}^{2})^{\sfrac{3}{2}}}.
\end{eqnarray}}\unskip

Here, $g$ is the acceleration of gravity, $T$ the surface tension
coefficient, $\rho _{1}$ and $\rho _{2}$ the densities. Herein, a
stable stratification is considered $(\mu < 1)$.
Equations~(\ref{eq3})  and~(\ref{eq4}) are, respectively, the kinematic
and the dynamic conditions.

On the solid boundaries, the problem is subjected to impermeability
conditions
{\begin{eqnarray}\label{eq8}
&\Phi _{z}^{(1)}=0 \quad Z=-{d}_{1}, \Seqnsplit
\label{eq9}
&{\Phi }_{z}^{(2)}= 0 \quad {Z}={d}_{2}, 
\end{eqnarray}}\unskip
with lower and upper boundaries at ${Z}=-{d}_{1}$, {and} ${Z}={d}_{2}$
respectively ($d_{1} > 0, d_2 >0$).

\section{Derivation of the algebraic equations}\label{sec3}

For irrotational interfacial waves, according to Whitham's
method~\cite{15}, we define the averaged Lagrangian over a rectangular
period, as 
{\begin{eqnarray}
\overline{{L}}&=&\frac{1}{4\uppi ^{2}}{{\int }}_{0}^{2\uppi }
{{\int }}_{0}^{2\uppi }\left\{\frac{1}{1-\mu }
{{\int }}_{-{d}_{1}}^{\eta }\left[-\varpi 
\frac{{\partial }\Phi _{1}}{{\partial }{X}}+
\frac{1}{2}\left(\overset{\rightarrow }
{{\nabla }}\Phi _{1}\right)^{2}\right]
\mathrm{d}Z\right\}\mathrm{d}X\mathrm{d}Y\nonumber\\
&& +\,\frac{1}{4\uppi ^{2}}{{\int }}_{0}^{2\uppi }
{{\int }}_{0}^{2\uppi }\left\{\frac{\mu }{1-\mu }
{{\int }}_{\eta }^{{d}_{2}}\left[-\varpi 
\frac{{\partial }\Phi _{2}}{{\partial }{X}}+
\frac{1}{2}\left(\overset{\rightarrow }{{\nabla }}
\Phi _{2}\right)^{2}\right]\mathrm{d}Z
+\frac{1}{2}\eta ^{2}+\kappa 
\left(\sqrt{1+\left(\overset{\rightarrow }
{{\nabla }}\eta \right)^{2}}-1\right)\right\}\mathrm{d}X\mathrm{d}Y
\label{eq10} \nonumber\\
\end{eqnarray}}\unskip
with
{\begin{equation}\label{eq11}
\overset{\rightarrow }{\nabla }=\left({p}\frac{\partial }
{\partial {X}},{q}\frac{\partial }{\partial {Y}},
\frac{\partial }{\partial {Z}}\right).
\end{equation}}\unskip

We seek expressions for $\Phi _{1}({X},{Y},{Z})$,  $\Phi
_{2}({X},{Y},{Z})$ and $\eta ({X},{Y})$  that are doubly periodic
functions of the transformed horizontal coordinates. A general
representation is 
{\begin{eqnarray}\label{eq12}
&\displaystyle\Phi _{1}({X},{Y},{Z})=\sum _{{mn}}{b}_{{mn}}
\frac{ch\alpha _{{mn}}({Z}+{d}_{1})}
{{ch\alpha }_{{mn}}{d}_{1}}\chi _{{mn}}({X},{Y}) \Seqnsplit
\label{eq13}
&\displaystyle\Phi _{2}({X},{Y},{Z})=\sum _{{mn}}{c}_{{mn}}
\frac{{ch\alpha }_{{mn}}({Z}-{d}_{2})}
{{ch\alpha }_{{mn}}{d}_{2}}\chi _{{mn}}({X},{Y}) \Seqnsplit
\label{eq14}
&\displaystyle\eta ({X,Y})=\sum_{m=0}^{\infty }
\sum_{n=0}^{\infty }\Delta _{{m0}}\Delta _{{n0}}{a}_{{mn}}
\cos mX \cos nY 
\end{eqnarray}}\unskip
where
{\begin{equation}\label{eq15}
\chi _{{mn}}({X,Y})=\sin mX \cos nY
\end{equation}}\unskip
is introduced for notational convenience to calculate the kinetic
energy and
{\begin{equation*}
\Delta _{j0}=1-\tfrac{1}{2}\delta _{j0}.
\end{equation*}}\unskip
Here ${\delta}_{j0}$ is the Kronecker symbol and
{\begin{equation}\label{eq16}
{\alpha }_{mn}^{2}=p^{2}m^{2}+q^{2}n^{2}.
\end{equation}}\unskip

Substituting expansions~(\ref{eq12})--(\ref{eq14})  into~(\ref{eq10}),
the averaged  Lagrangian $\overline{L}$ is obtained as a function of
the  unknown coefficients ${a_{mn}}$, $b_{mn}$ and ${c}_{{mn}}$. Some
details of the calculation of $\overline{L}$  are given in the
Appendix~\ref{secA}. The resulting averaged Lagrangian is expressed in
the following way
{\begin{eqnarray}\label{eq17}
\overline{{L}}&=&\frac{1}{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{b}_{{mn}}{b}_{{kl}}}{{ch}
(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times \frac{1}
{(\alpha _{{mn}}+\alpha _{{kl}})}
[\mbox{}_{1}^{}{\Omega }{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}\neq 
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
[\mbox{}_{1}^{}{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}{{ch}^{2}
(\alpha _{{mn}}{d}_{1})}\times 
[\mbox{}_{1}^{} {\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times \frac{1}
{(\alpha _{{mn}}+\alpha _{{kl}})}
[\mbox{}_{2}^{}{\Omega }{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}
\neq \alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
[\mbox{}_{2}^{}{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}^{2}(\alpha _{{mn}}{d}_{2})}\times 
[\mbox{}_{2}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\nonumber\\
&&+\,\frac{1}{8}\sum_{{m}=0}^{\infty }\sum_{{n}=0}^{\infty }
\Delta _{{m}0}\Delta _{{n}0}{{a}}_{{mn}}^{2}-
\frac{1}{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{m}\frac{{b}_{{mn}}}{\alpha _{{mn}}{ch}
(\alpha _{{mn}}{d}_{1})}
\mbox{}_{1}^{}
{\omega }{_{mn}^{mn}}\nonumber\\
&&+\,\frac{\mu }{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{m}\frac{{c}_{{mn}}}{\alpha _{{mn}}
{ch}(\alpha _{{mn}}{d}_{2})}
\mbox{}_{2}^{}
{\omega }{_{mn}^{mn}}+\frac{\kappa }{4}
({R}_{00}-4).
\end{eqnarray}}\unskip

The quantities $[{\Theta }_{{mnkl}}^{+}],[{\Theta
}_{{mnkl}}^{-}],[\mbox{}_{i}^{}{\Omega
}{_{mnkl}^{+}}],[\mbox{}_{i}^{}{\Omega
}{_{mnkl}^{-}}],[\mbox{}_{i}^{}{\Omega }{_{mn}^{}}]$ and
$[\mbox{}_{i}^{}{\Lambda }{_{mn}^{}}]({i}=1,2)$, are defined in the
Appendix~\ref{secA}. 


According to Whitham's theory, the variational equations to be solved
are
{\begin{eqnarray}\label{eq18}
&\displaystyle\frac{\partial \overline{L}}{\partial a_{mn}}=0 \Seqnsplit
\label{eq19} 
&\displaystyle\frac{\partial \overline{L}}{\partial b_{mn}}=0.
\end{eqnarray}}\unskip

Applying Equations~(\ref{eq18}) and~(\ref{eq19}) to the averaged
Lagrangian~(\ref{eq17}) yields the following set of nonlinear
algebraic equations 
{\begin{eqnarray}
\frac{\partial \overline{{L}}}{\partial {b}_{{rs}}}
&=&\frac{1}{32(1-\mu)}\sum_{{m}=1}^{\infty }\sum_{{n}=0}^{\infty }
\frac{{b}_{{mn}}}{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{rs}}{d}_{1})}\times \frac{1}{(\alpha _{{mn}}+\alpha _{{rs}})}
[\mbox{}_{1}^{}{\Omega }{_{mnrs}^{+}}]^{ {t}} 
[{\Theta }_{{mnsr}}^{+}]\nonumber\\
&&+\,\frac{1}{32(1-\mu)}\underset{\alpha _{{mn}}\neq 
\alpha _{{rs}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }}}\frac{{b}_{{mn}}}{{ch}
(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{rs}}{d}_{1})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{rs}})}
[\mbox{}_{1}^{}{\Omega }{_{mnrs}^{-}}] ^{{t} }
[{\Theta }_{{mnrs}}^{-}]\nonumber\\
&&+\,\frac{1}{32(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{rs}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }}}\frac{{b}_{{mn}}}{{ch}^{2}
(\alpha _{{rs}}{d}_{1})}\times 
[\mbox{}_{1}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{1}{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\frac{{r}}{\alpha _{{rs}}{ch}
(\alpha _{{rs}}{d}_{1})}\mbox{}_{1}^{}
{\omega }{_{rs}^{rs}}=0\label{eq20}\Seqnsplit
\frac{\partial \overline{{L}}}{\partial {c}_{{rs}}}&=&
\frac{1}{32(1-\mu)}\sum_{{m}=1}^{\infty }\sum_{{n}=0}^{\infty }
\frac{{c}_{{mn}}}{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{rs}}{d}_{2})}\times \frac{1}
{(\alpha _{{mn}}+\alpha _{{rs}})}
[\mbox{}_{2}^{}{\Omega }{_{mnrs}^{+}}]^{ {t}} 
[{\Theta }_{{mnsr}}^{+}]\nonumber\\
&&+\,\frac{1}{32(1-\mu)}\underset{\alpha _{{mn}}\neq 
\alpha _{{rs}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }}}\frac{{c}_{{mn}}}{{ch}
(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{rs}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{rs}})}
[\mbox{}_{2}^{}{\Omega }{_{mnrs}^{-}}] ^{{t} }
[{\Theta }_{{mnrs}}^{-}]\nonumber\\
&&+\,\frac{1}{32(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{rs}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }}}\frac{{c}_{{mn}}}{{ch}^{2}
(\alpha _{{rs}}{d}_{2})}\times 
[\mbox{}_{2}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{1}{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\frac{{r}}{\alpha _{{rs}}{ch}
(\alpha _{{rs}}{d}_{2})}\mbox{}_{2}^{}
{\omega }{_{rs}^{rs}}=0\label{eq21}\Seqnsplit
\frac{\partial \overline{{L}}}{\partial {a}_{{rs}}}&=&
\frac{1}{64(1-\mu)}\sum_{{m}=1}^{\infty }\sum_{{n}=0}^{\infty }
\sum_{{k}=1}^{\infty }\sum_{{l}=0}^{\infty }
\frac{{b}_{{mn}}{b}_{{kl}}}{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times \frac{1}
{(\alpha _{{mn}}+\alpha _{{kl}})}
\frac{\partial }{\partial {a}_{{rs}}}[\mbox{}_{1}^{}
{\overline{\Omega }}{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\vspace*{3pt}\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}\neq 
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}{{ch}
(\alpha _{{mn}}{d}_{1}){sh}
(\alpha _{{kl}}{d}_{1})}\times \frac{1}
{(\alpha _{{mn}}-\alpha _{{kl}})}
\frac{\partial }{\partial {a}_{{rs}}}[\mbox{}_{1}^{}
{\Omega }{_{mnkl}^{-}}] ^{{t} }[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}=\alpha _{{kl}}}
{\underbrace{\sum_{{m}=1}^{\infty }\sum_{{n}=0}^{\infty }
\sum_{{k}=1}^{\infty }\sum_{{l}=0}^{\infty }}}
\frac{{b}_{{mn}}{b}_{{kl}}}{{ch}^{2}(\alpha _{{mn}}{d}_{1})}
\frac{\partial }{\partial {a}_{{rs}}}[\mbox{}_{1}^{}
{\Lambda }{_{mn}^{}}]^{ {t}} [{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{c}_{{mn}}{c}_{{kl}}}{{ch}
(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
\frac{\partial }{\partial {a}_{{rs}}}[\mbox{}_{2}^{}
{\Omega }{_{mnkl}^{+}}]^{ {t}} [{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}\neq 
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times \frac{1}
{(\alpha _{{mn}}-\alpha _{{kl}})}
\frac{\partial }{\partial {a}_{{rs}}}[\mbox{}_{2}^{}
{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}^{2}(\alpha _{{mn}}{d}_{2})}\times 
\frac{\partial }{\partial {a}_{{rs}}}
[\mbox{}_{2}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{1}{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{m}\frac{{b}_{{mn}}}{\alpha _{{mn}}{ch}
(\alpha _{{mn}}{d}_{1})}\frac{\partial }
{\partial {a}_{{rs}}}(\mbox{}_{1}^{}{\omega }{_{mn}^{mn}}
)\nonumber\\
&&+\,\frac{\mu }{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{m}\frac{{c}_{{mn}}}{\alpha _{{mn}}
{ch}(\alpha _{{mn}}{d}_{2})}\frac{\partial }
{\partial {a}_{{rs}}}(\mbox{}_{2}^{}
{\omega }{_{mn}^{mn}})\nonumber\\
&&+\,\frac{\kappa }{16}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{a}_{{mn}}|\overline{{T}}_{{mnrs}}| ^{{t}}
| \vartheta _{\text{mnrs}}| +\frac{1}{4}\Delta _{{r}0}
\Delta _{{s}0}{a}_{{rs}}=0.
\label{eq22}
\end{eqnarray}}\unskip

\section{Numerical procedure}\label{sec4}

The series in~(\ref{eq12})--(\ref{eq14}) are truncated at N and the
order N is chosen for truncating the other expansions and the series in
Equations~(\ref{eq20})--(\ref{eq22}). The resulting system 
is to be solved for the unknown 
$a_{mn}, b _{mn}, c _{mn}$ and ${\varpi}$ 
by fixing the values of ${\theta}$, ${\mu}$, 
$d_{1}, d_{2}$ and the steepness which is defined as 
{\begin{equation}\label{eq23}
\varepsilon =\eta (0,0)-\eta (\pi,0)=
\sum_{{m}=0}^{{N}} \sum_{{n}=0}^{{N}}{a}_{2{m}+1,2{n}+1}.
\end{equation}}\unskip
Of course, this definition is different from the peak-to-through wave
height, since the maximum and minimum points are not always at
$(x,y)=(0,0)$ and $(\pi,0)$, respectively. 

Note that $a_{00}=0$, the mean surface level being at $Z=0$ and, owing to
the triangular symmetry of a short-crested wave,
{\begin{equation}\label{eq24}
a_{{mm}}={b}_{{mn}}={c}_{{mm}}=0\quad  \text{when } m+n \text{ is odd}.
\end{equation}}\unskip

The system of equations~(\ref{eq20})--(\ref{eq23}) is solved by Newton's
method. In many practical computations it has been found that less than
five iterations are sufficient to satisfy the system of equations with
an error less than $10^{-14}$. Of course, the required number of
iterations depends on the configuration studied and on the desired
accuracy. For instance, this is the case when we have a slow decay of
the Fourier coefficients. To get convergent series it is necessary to
increase the number of coefficients, which causes higher computational
coasts. Decreasing the order of the truncated series often results in
increasing the number of iterations.

The initial estimates for the Newton iteration are given for small
amplitudes by analytical solutions calculated up to the  ``fifth-order''
using a perturbation procedure and the computer algebra system MAPLE
for symbolic computation. 

The Fourier coefficients introduced for the purposes of the present
method are computed using the fast Fourier transform (FFT) algorithm.
The computational domain is rectangular, periodic in both directions
and discretized with constant steps $\Delta X=\Delta Y=2\uppi /M$. $M$
and $N$ may be chosen with a weak constraint, i.e. $M > 2N$, in order to
avoid the aliasing phenomenon. 


To check the accuracy of the present method, it was necessary to
compare our results with those obtained, in the case of gravity
interfacial waves, by Allalou \etal~\cite{12}. A comparison between
computed values of the phase velocity is presented in Table~\ref{tab1}. Results
using the method just described are in excellent agreement to those of
Allalou \etal~computed via a perturbation expansion up to the 27th
order in wave steepness.

\begin{table}
\caption{\label{tab1}Comparison between computed values of the phase
velocity with the present method and those given by 
Allalou \etal~\cite{12}
for $\varepsilon =0.01,{d}_{1}=0.8,{d}_{2}=1,\theta
=50{\text{\textdegree}},\mu =0.1$ and $0.5$}
\begin{tabular}{ccc}
\thead
${\mu}$ & Present method ($N=23$) & Perturbation method ($N=27$) \\
\endthead
0.1 & 0.7443973837\0 &  0.7443973838\0  \\
0.5 &   0.48129554183  &    0.48129554186 
\botline
\end{tabular}
\end{table}

\section{Results and discussion}\label{sec5}

To illustrate the performance of the method we shall present examples
of fully three-dimensional calculations. An interesting variety of wave
patterns emerges, with a distinct change in character as a result of
different combinations of wave parameters.

\subsection{Resonant cases}\label{sec51}

To find the equation which relates the parameters of short crested
waves for which harmonic resonance occurs with respect to the ($m,
n$)th harmonic, we solve the system~(\ref{eq2})--(\ref{eq4}) by using
the classical perturbation method. The perturbation parameter we use is
the steepness ${\varepsilon}$. We expand the quantities $\Phi
_{{1}},\Phi _{{2}},\eta $ {and ${\varpi}$} as power series
in the parameter ${\varepsilon}$.
{\begin{equation}\label{eq25}
\Phi _{{1}}=\sum _{{r}={1}}{\Phi }_{{1}}^{({r})}
\varepsilon ^{{r}}\quad\Phi _{{2}}=\sum _{{r}={1}}
{\Phi }_{{2}}^{({r})}\varepsilon ^{{r}}\quad
\eta =\sum _{{r}={1}}\eta ^{({r})}\varepsilon ^{{r}}\quad
\varpi =\sum _{{r}={0}}\varpi _{{r}}\varepsilon ^{{r}}.
\end{equation}}\unskip

The dimensionless velocity potentials at the free surface are expressed
in terms of the Taylor expansions at $Z=0$ instead of ${Z}=\eta
({X,Y})$:
{\begin{equation}\label{eq26}
\Phi _{{i}}({X,Y,Z}=\eta ({X,Y}))=\Phi _{{i}}({X,Y,0})+
\left.\frac{\partial \Phi _{{i}}}{\partial {Z}}\right| _{{Z}={0}}
\eta +\left.\frac{\partial ^{{2}}\Phi _{{i}}}{\partial {Z}^{{2}}}
\right| _{{Z}={0}}\frac{\eta ^{{2}}}{{2}!}+\left.\frac{\partial ^{{3}}
\Phi _{{i}}}{\partial {Z}^{{3}}}\right| _{{Z}={0}}
\frac{\eta ^{{3}}}{{3}!}+\cdots .
\end{equation}}\unskip
As did Kimmoun \etal~\cite{5}, ${\xi}(X,Y)$ is written in this form
{\begin{equation}\label{eq27}
\xi ({X},{Y})=\frac{{N}({X},{Y})}
{(1+{D}({X},{Y}))^{\sfrac{3}{2}}}
\end{equation}}\unskip
with
{\begin{equation*}
{N}({X},{Y})={p}^{2}\eta _{{XX}}(1+{q}^{2}{\eta }_{{Y}}^{2})+
{q}^{2}\eta _{{YY}}(1+{p}^{2}{\eta }_{{X}}^{2})-
2{p}^{2}{q}^{2}\eta _{{XY}}\eta _{{X}}\eta _{{Y}}
\end{equation*}}\unskip
and
{\begin{equation*}
{D}({X},{Y})={p}^{2}{\eta }_{{X}}^{2}+{q}^{2}{\eta }_{{Y}}^{2}
\end{equation*}}\unskip
${\xi}(X,Y)$ is then developed in terms of a series expansion:
{\begin{equation}\label{eq28}
\xi ({X},{Y})={N}({X},{Y})\left[1-\tfrac{3}{2}{D}+
\tfrac{15}{8}{D}^{2}-\tfrac{105}{48}{D}^{3}+{O}({D}^{5})\right].
\end{equation}}\unskip
After substitution into Equations~(\ref{eq2})--(\ref{eq4}) and
gathering the ${\varepsilon}$ of the same power, we get a new set of
equations that was processed using the symbolic software package MAPLE.
We proceed by seeking solutions of the form
{{\begin{equation}\label{eq29}
\begin{cases} 
\displaystyle{\Phi }_{{1}}^{({r})}=\sum \sum {{b}}_{{mn}}^{({r})}
\sin ({mX})\cos ({nY})\frac{\cosh 
(\alpha _{{mn}}({Z}+{d}_{{1}}))}
{\cosh (\alpha _{{mn}}{d}_{{1}})}\vspace*{5pt}\\
\displaystyle{\Phi }_{{2}}^{({r})}=\sum \sum {{c}}_{{mn}}^{({r})}\sin 
({mX})\cos ({nY})\frac{\cosh (\alpha _{{mn}}
({Z}-{d}_{{2}}))}{\cosh (\alpha _{{mn}}{d}_{{2}})}\vspace*{5pt}\\
\displaystyle\eta ^{({r})}=\sum \sum {{a}}_{{mn}}^{({r})}\cos ({mX})\cos ({nY}).
\end{cases} 
\end{equation}}\unskip
To determine the resonance condition, one can rewrite the kinematic and
dynamical equations, at order ($r$), as follows
{\begin{equation}
\label{eq30}
\begin{cases} 
\displaystyle-{m\varpi }_{{0}}{{a}}_{{mn}}^{({r})}+
\alpha _{{mn}}{F}_{{mn}}{{b}}_{{mn}}^{({r})}={{A}}_{{mn}}^{({r}-{1})}\vspace*{5pt}\\
\displaystyle-{m\varpi }_{{0}}{{a}}_{{mn}}^{({r})}-
\alpha _{{mn}}{G}_{{mn}}{{c}}_{{mn}}^{({r})}={{B}}_{{mn}}^{({r}-{1})}\vspace*{5pt}\\
\displaystyle({\kappa \alpha }_{{mn}}^{{2}}+{1}-\mu)
{{a}}_{{mn}}^{({r})}+{m\varpi }_{{0}}
({{\mu c}}_{{mn}}^{({r})}-{{b}}_{{mn}}^{({r})})=
{{C}}_{{mn}}^{({r}-{1})}
\end{cases} 
\end{equation}}\unskip
where $A_{{mn}}^{({r}-{1})}$, ${{B}}_{{mn}}^{({r}-{1})}$ {and}
${{C}}_{{mn}}^{({r}-{1})}$  are terms depending only on coefficients
of order ($r-1$), 
{\begin{equation}\label{eq31}
{F}_{{mn}}=\tanh (\alpha _{{mn}}{d}_{1}),\quad
{G}_{{mn}}=\tanh (\alpha _{{mn}}{d}_{2})
\quad\text{and}\quad {\varpi }_{0}^{2}=\frac{(1-\mu)
(1+\kappa)\tanh ({d}_{1})\tanh ({d}_{2})}{\mu \tanh({d}_{2})+
\tanh ({d}_{1})}.
\end{equation}}\unskip
The above system of equations has solutions of the form
{\begin{equation}
\label{eq32}
\begin{cases} 
\displaystyle{{a}}_{{mn}}^{({r})}=-\frac{\alpha _{{mn}}{G}_{{mn}}{F}_{{mn}}
{{C}}_{{mn}}^{({r}-{1})}+{m\varpi }_{{0}}
\left({G}_{{mn}}{{A}}_{{mn}}^{({r}-{1})}+{\mu F}_{{mn}}{
{B}}_{{mn}}^{({r}-{1})}\right)}{-\alpha _{{mn}}
({\kappa \alpha }_{{mn}}^{{2}}+{1}-\mu)
{F}_{{mn}}{G}_{{mn}}+{m}^{{2}}{\varpi }_{{0}}^{{2}}
({G}_{{mn}}+{\mu F}_{{mn}})}\vspace*{5pt}\\
\displaystyle{{b}}_{{mn}}^{({r})}=\frac{{m}^{{2}}{\varpi }_{{0}}^{{2}}
\mu \left({{A}}_{{mn}}^{({r}-{1})}-{
{B}}_{{mn}}^{({r}-{1})}\right)-\alpha _{{mn}}{m\varpi }_{{0}}
{G}_{{mn}}{{C}}_{{mn}}^{({r}-{1})}-\alpha _{{mn}}{G}_{{mn}}
{{A}}_{{mn}}^{({r}-{1})}({\kappa \alpha }_{{mn}}^{{2}}+
{1}-\mu )}{-\alpha _{{mn}}({\kappa \alpha }_{{mn}}^{{2}}
+{1}-\mu){F}_{{mn}}{G}_{{mn}}+{m}^{{2}}{\varpi }_{{0}}^{{2}}
({G}_{{mn}}+{\mu F}_{{mn}})}\vspace*{5pt}\\
\displaystyle{{c}}_{{mn}}^{({r})}=\frac{{m}^{{2}}{\varpi }_{{0}}^{{2}}
\left({{A}}_{{mn}}^{({r}-{1})}-{{B}}_{{mn}}^{({r}-{1})}\right)+
\alpha _{{mn}}{m\varpi }_{{0}}{F}_{{mn}}{{C}}_{{mn}}^{({r}-{1})}+
\alpha _{{mn}}{F}_{{mn}}{{B}}_{{mn}}^{({r}-{1})}
({\kappa \alpha }_{{mn}}^{{2}}+{1}-\mu)}
{-\alpha _{{mn}}({\kappa \alpha }_{{mn}}^{{2}}+{1}-
\mu ){F}_{{mn}}{G}_{{mn}}+{m}^{{2}}{\varpi }_{{0}}^{{2}}
({G}_{{mn}}+{\mu F}_{{mn}})}.
\end{cases} 
\vspace*{5pt}
\end{equation}}\unskip
}

These coefficients have zero divisors when the parameters of the wave
satisfy the relation   \vspace*{5pt}
{\begin{equation}\label{eq33}
-\alpha _{{mn}}({\kappa \alpha }_{{mn}}^{{2}}+{1}-
\mu)\tanh (\alpha _{{mn}}{d}_{1})
\tanh (\alpha _{{mn}}{d}_{2})+{m}^{{2}}
{\varpi }_{{0}}^{{2}}\left(\tanh (\alpha _{{mn}}
{d}_{2})+\mu \tanh(\alpha _{{mn}}{d}_{1})\right)=0.\vspace*{5pt}
\end{equation}}\unskip
These zero divisors are now recognized to be indicative of the
occurrence of harmonic resonances. Physically, the fundamental ($1, 1$)
excites the harmonic ($m, n$) which travels at the same phase speed. Over
a long time scale the resonance allows continual redistribution of
energy between the fundamental and the resonant harmonic. 

Taking advantage of the ability of the present method to provide the
properties of the fully resonant waves, it is natural to examine these
resonances. Taking ${m}={n}=2$, ${d}_{1}={d}_{2}=0.5,\theta
=48{\text{\textdegree}}$ and $\mu =0.07$,  Equation~(\ref{eq33}) yields
$\kappa =0.07587182693956657\ldots.$ As known, such a resonance
gives rise to a multiplicity of solutions. In our case, we obtain three
solution branches; two of them being connected through a turning point
which occurs at about ${\varepsilon}=0.009$, as shown in
Figure~\ref{fig2}. The latter exhibits the multiple solution structure
of the wave speed as a function of steepness. The starting point of
branch~(\ref{eq1}) was generated from a non resonant solution computed
with $\kappa =0.08$ varying progressively ${\kappa}$ from this
value to $0.07587182693956657$ by step not exceeding few percents of
the fixed value. Figure~\ref{fig3} shows the shape of this family of
{\vadjust{\pagebreak}}short-crested waves. The lines of the intersection with the planes
${X}=2\uppi$ or $Y=0$ (remind that the patterns are 
$2\uppi$-periodic in the both directions) have the same shape as the profile
of Wilton type 1 ripples~\cite{16}: they exhibit secondary troughs in
the crests. The starting point of branch~(\ref{eq2}) was obtained, 
for small values of the wave steepness, by taking 
$\kappa=0.07587182693956657$.  Figure~\ref{fig4} 
shows that the lines of the
intersection of their shape with the planes ${X}=2\uppi$ 
or $Y=0$ are similar to the profile of the so-called Wilton
type 2 ripples, that is to say, secondary crests of a substantial
height appear in the troughs. Beyond the turning point, we have a
mixed-type solution which combines the latter two patterns.
Figure~\ref{fig5} shows the shape of these waves which exhibits
secondary troughs in the crests and similarly secondary crests in the
troughs. 

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}Variation of wave speed with steepness for
three-dimensional Wilton ripples.}
\end{figure}

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}Three-dimensional Wilton ripples of the type 1
for ${d}_{1}={d}_{2}=0.5$, $\theta =48{\text{\textdegree}}$, $\mu =0.07$
and $\kappa =0.07587182693956657$.}
{\vspace*{4pt}}
\end{figure}

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig4}Three-dimensional Wilton ripples of type 2 for
${d}_{1}={d}_{2}=0.5$, $\theta =48{\text{\textdegree}}$, $\mu =0.07$ and
$\kappa =0.07587182693956657$.}
\end{figure}

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig5}Three-dimensional Wilton ripples of type 3 for
${d}_{1}={d}_{2}=0.5$, $\theta =48{\text{\textdegree}}$, $\mu =0.07$ and
$\kappa =0.07587182693956657$.}
\end{figure}

\subsection{Energies}\label{sec52}

Once the coefficients $a_{mn}$, $b_{mn}$ and $c_{mn}$  determined, it
is straightforward to compute kinetic and potential energies by using
the expressions~(\ref{eqA9}),~(\ref{eqA10}) and~(\ref{eqA13}) given in
the Appendix~\ref{secA}.  As an illustration we present the results
obtained for the resonant cases discussed above. The mean kinetic-~and
potential-energy densities for the wave referred to as type 1 are
plotted in Figure~\ref{fig6}. Those of the cases referred to as type 2
and type 3 are plotted in Figure~\ref{fig7}. In Figure~\ref{fig8}, the
energy densities of the wave of type 2 correspond to
$\varepsilon<0.009$ and those of the wave of type 3 to $\varepsilon
>0.009$. For all the cases presented in these two figures, the
kinetic-energy density is larger than the two potential-energy
densities. There are configurations for which potential-energy density,
${V}_{{g}}$, is greater than the kinetic-energy, as shown in
Figure~\ref{fig8}. However, we did not find any cases where the
capillary-potential-energy density, ${V}_{\tau}$, is greater than
the other two energy densities.

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig6}Variation of mean kinetic energy density ($KE$),
mean gravity-potential-energy density $(V_g)$ and mean
capillary-potential-energy density with steepness for
three-dimensional Wilton ripples of type 1.}
\end{figure}

\begin{figure}
\includegraphics{fig07}
\caption{\label{fig7}Variation of mean kinetic energy density ($KE$),
mean gravity-potential-energy density ($V_g$) and mean
capillary-potential-energy density with steepness for three-dimensional
Wilton ripples of the type 2 and of the type 3.}
\end{figure}

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig8}Variation of mean kinetic energy density ($KE$) and
mean gravity-potential-energy density $(V_{g})$ with the density
ratio ${\mu}$ for ${d}_{1}={d}_{2}=0.1$, $\theta =50{\text{\textdegree}}$,
$\varepsilon =0.02$  and $\kappa =0$.}
\end{figure}

\subsection{Square waves}\label{sec53}

When one layer is shallow and the other is deep, the wave profile is
generally more sharply peaked on the side pointing into the deeper
layer. When the lower fluid has a very small depth and the upper fluid
is deep, it will be possible to obtain strongly localized shapes
similar to grid waves that form square-shaped patterns like those shown
in Figure~\ref{fig9}. A similar surface wave pattern has been
photographed by Terry Toedte-Meier off the coast of Oregon 
(cf.~Segur~\cite{17}). Solutions in the form of grid waves can be treated as
a composition of two plane solitary waves propagating at right angles
to each other. Each wave train exhibits flat-trough, sharp-crest
pattern, but successive wave crests in the same train are so far apart
that each crest behaves nearly like a solitary wave. At the
intersections of the crests appear peaks of greater elevation. The
flatness of the trough is due to effects from the bottom and may
indicate that these solutions are more nonlinear. If, on the contrary,
it is the upper layer that is shallow and the lower layer is deep, we
obtain similar  square-shaped patterns, but their convex sides are
oriented upward as shown in Figure~\ref{fig10}. These two families of
waves are obtained more easily with very low values of ${\mu}$ and
${\kappa}$. Their chapes changes as we increase these two parameters.

\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig09}
\vspace*{-.5pc}
\caption{\label{fig9}Graph in perspective of square waves with their
convex sides oriented downward for ${d}_{1}=0.034, {d}_{2}=0.1$,
$\theta =1{\text{\textdegree}}$, $\varepsilon =0.02,\mu =0.1$ and $\kappa
=0$. The vertical scale is enlarged to emphasize the shape.}
{\vspace*{-3pt}}
\end{figure}

\begin{figure}
\includegraphics{fig10}
\vspace*{-.5pc}
\caption{\label{fig10}Graph in perspective of square waves with their
convex sides oriented upward for ${d}_{1}=5$, ${d}_{2}=0.02$, $\theta
=1{\text{\textdegree}}$, $\varepsilon =0.02$, $\mu =0.1$ and $\kappa =0$.
The vertical scale is enlarged to emphasize the shape.}
{\vspace*{-3pt}}
\end{figure}

\subsection{Hexagonal waves}\label{sec54}

For low values of ${\theta}$, our method can generate hexagonal
traveling wave solutions. There exist several types of these waves.
When the lower fluid has a very small depth and the upper fluid is
deep, we can generate honeycomb shapes like those shown in
Figure~\ref{fig11}. Each cell have a large flat trough which is
surrounded by a narrow six sided ridge. The pattern is elongated in the
direction perpendicular to the phase velocity. The sides of the ridge
which are aligned with the $Y$-direction are relatively flat over a
large fraction of their length and are higher than the other four sides
which complete the hexagonal pattern. If it is the upper layer that is
shallow and the lower layer is deep, we obtain similar hexagonal
patterns but each of them have a large flat crest which is surrounded
by a six sided depression boundary (see Figure~\ref{fig12}).
Figure~\ref{fig13} presents a third notable case example which exhibit
mixed structures, that is to say they have well-pronounced diamond
crests while the troughs have hexagonal patterns. Hexagonal patterns
have been observed experimentally by Hammack \etal~\cite{18}.

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig11}Graph in perspective of hexagonal waves with
their convex sides oriented downward for ${d}_{1}=0.034$, ${d}_{2}=0.1$,
$\theta =5{\text{\textdegree}}$, $\varepsilon =0.02$, $\mu =0.1$ and
$\kappa =0$. The vertical scale is enlarged to emphasize the
shape.}
\end{figure}

\begin{figure}
{\vspace*{6pt}}
\includegraphics{fig12}
{\vspace*{6pt}}
\caption{\label{fig12}Graph in perspective of hexagonal waves with
their convex sides oriented upward for ${d}_{1}=3$, ${d}_{2}=0.02$,
$\theta =50{\text{\textdegree}}$, $\varepsilon =0.02$, $\mu =0.1$ and $\kappa
=0$. The vertical scale is enlarged to emphasize the shape.}
{\vspace*{6pt}}
\end{figure}

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig13}Graph in perspective of mixed waves.}
\end{figure}

\section{Conclusion}\label{sec6}

A method has been developed for the solution of the problem involving
interfacial capillary--gravity short-crested waves propagating at the
interface between two inviscid immiscible fluids of different densities
lying between two rigid boundaries. A large number of interesting
configurations can be computed using our approach. Several results have
been presented to illustrate the effectiveness of our approach. The
shapes of these waves depend on different parameters of the problem:
the surface tension coefficient, the fluid layer depth, density
ratios and the angle between the direction of propagation of the
incident wave and the normal to the wall. The influence of these
parameters will be further investigated in the future. 

\section*{Conflicts of interest}

Authors have no conflict of interest to declare.

\section*{Acknowledgements}

The authors would like to thank Professor A. Ouibrahim for his invaluable
suggestions and for his helpful remarks to the text.

\back{}


\appendix{}

\def\appendixlabel{Appendix A}

\setcounter{equation}{0}
\renewcommand\theequation{A\arabic{equation}}


\section{Calculation of the averaged Lagrangian $\overline{L}$}\label{secA}


This appendix describes some details of the calculation of the averaged
Lagrangian $\overline{L}$.

\subsection{The term of kinetic energy}\label{secA1}

In the averaged Lagrangian function~(\ref{eq10}) given in the body of
the paper, the contribution of the kinetic energy is
{\begin{equation}\label{eqA1}
\mathrm{KE}= \underset{\mathrm{KE}_{{1}}}{\underbrace{\frac{{1}}{{1}-\mu }
\frac{{1}}{{8\uppi }^{{2}}}\int_{{0}}^{{2\uppi }} 
\int_{{0}}^{{2\uppi }} \int_{-{d}_{{1}}}^{\eta }
\left(\overset{\rightarrow }{\nabla }\Phi _{{1}}\right)^{{2}}
\mathrm{d}X  \mathrm{d}Y \mathrm{d}Z}}+\underset{\mathrm{KE}_{{2}}}
{\underbrace{\frac{\mu }{{1}-\mu }\frac{{1}}
{{8\uppi }^{{2}}}\int_{{0}}^{{2\uppi }} 
\int_{{0}}^{{2\uppi }} \int_{\eta }^{{d}_{{2}}}
\left(\overset{\rightarrow }{\nabla }\Phi _{{2}}\right)^{{2}}
\mathrm{d}X  \mathrm{d}Y \mathrm{d}Z}}.
\end{equation}}\unskip

The first term on the right-hand side of Equation~(\ref{eqA1}), 
labelled $\mathrm{KE}_{{1}}$, represents the kinetic energy in the lower
fluid of the wave. Substituting expansion~(\ref{eq12}) into 
$\mathrm{KE}_{{1}}$ leads to
{\begin{eqnarray}\label{eqA2}
\mathrm{KE}_{{1}}&=&\frac{1}{1-\mu }\frac{1}{8\uppi ^{2}}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\left(\frac{{b}_{{mn}}{b}_{{kl}}}{{ch}
(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\right) \nonumber\\ 
&&\times\,\int_{0}^{2\uppi } \int_{0}^{2\uppi } \left\{\int_{-{d}_{1}}^{\eta }
{ch}(\alpha _{{mn}}({Z}+{d}_{1}))
{ch}(\alpha _{{kl}}({Z}+{d}_{1}))
\mathrm{d}Z\,\overset{\rightarrow }{\nabla }(\chi _{{mn}})
\cdot \overset{\rightarrow }{\nabla }
(\chi _{{kl}})\right.\nonumber\\ 
&& +\,\alpha _{{mn}}\alpha _{{kl}}\left.\int_{-{d}_{1}}^{\eta }
{sh}(\alpha _{{mn}}({Z}+{d}_{1}))
{sh}(\alpha _{{kl}}({Z}+{d}_{1}))
\mathrm{d}Z\,\chi_{{mn}}\chi _{{kl}}\right\}\mathrm{d}X \mathrm{d}Y.
\end{eqnarray}}\unskip

Using the properties of the products of hyperbolic functions and
proceeding in a similar way as in Debiane and Kharif~\cite{11}, we
obtain
{\begin{eqnarray}\label{eqA3}
&&\mathrm{KE}_{{1}}=\frac{1}{1-\mu }\frac{1}{64}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})} \nonumber\\
&&\,\,\,\times\, \frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
\left\{\!\begin{array}{@{}l@{}}
\left[\mbox{}_{1}^{}{{P}}{_{{k}-{m},{l}-{n}}^{mnkl}}({p}^{2}{mk}+
{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})+\mbox{}_{1}^{}
{{P}}{_{{k}+{m},{l}-{n}}^{mnkl}}({p}^{2}{mk}-{q}^{2}{nl}-
\alpha _{{mn}}\alpha _{{kl}})\right.\\
\left.+\,\mbox{}_{1}^{}{{P}}{_{{k}-{m},{l}+{n}}^{mnkl}}
({p}^{2}{mk}-{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})+
\mbox{}_{1}^{}{{P}}{_{{k}+{m},{l}+{n}}^{mnkl}}
({p}^{2}{mk}+{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})\right]\end{array}\!\right\}
\nonumber\\
&&\,\,\,+\,\frac{1}{1-\mu }\frac{1}{64}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})} \nonumber\\
&&\,\,\,\times\, \frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
\left\{\!\begin{array}{@{}l@{}}
\left[\mbox{}_{1}^{}{{M}}{_{{k}-{m},{l}-{n}}^{mnkl}}
({p}^{2}{mk}+{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})+
\mbox{}_{1}^{}{{M}}{_{{k}+{m},{l}-{n}}^{mnkl}}({p}^{2}{mk}-
{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})\right.\\
\left.+\,\mbox{}_{1}^{}{{M}}{_{{k}-{m},{l}+{n}}^{mnkl}}
({p}^{2}{mk}-{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})+
\mbox{}_{1}^{}{{M}}{_{{k}+{m},{l}+{n}}^{mnkl}}({p}^{2}{mk}+
{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})\right]\end{array}\!\right\}\nonumber\\
\end{eqnarray}}\unskip
$\mbox{}_{1}^{}{{P}}{_{rs}^{mnkl}}$ and
$\mbox{}_{1}^{}{{M}}{_{rs}^{mnkl}}$ being Fourier coefficients
defined by
{\begin{eqnarray}
{sh}(\alpha _{{mn}}+\alpha _{{kl}})
(\eta +{d}_{1})&=&
\sum_{{u}=0}^{{\infty }}\sum_{{v}=0}^{{\infty }}
\Delta _{{u}0}\Delta _{{v}0}\,\mbox{}_{1}^{}
{{P}}{_{uv}^{mnkl}} \cos uX \cos vY\label{eqA4}\Seqnsplit
{sh}(\alpha _{{mn}}+\alpha _{{kl}})
(\eta +{d}_{1})&=&\sum_{{u}=0}^{{\infty }}\sum_{{v}=0}^{{\infty }}
\Delta _{{u}0}\Delta _{{v}0}\mbox{}_{1}^{}{{M}}{_{uv}^{mnkl}} 
\cos uX \cos vY.\label{eqA5}
\end{eqnarray}}\unskip

The second term on the right-hand side of Equation~(\ref{eqA1}),
labelled $\mathrm{KE}_{2}$, represents the kinetic energy in the upper fluid
of the wave. This integral is handled in exactly the same way as for
$\mathrm{KE}_{{1}}$. The result is
{\begin{eqnarray}
&&\mathrm{KE}_{2}=-\frac{\mu }{1-\mu }\frac{1}{64}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{c}_{{mn}}{c}_{{kl}}}{{ch}
(\alpha _{{mn}}{d}_{2}){ch}(\alpha _{{kl}}{d}_{2})}\nonumber\\
&&\,\,\,\,\times\, \frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
\left\{\!\begin{array}{@{}l@{}}\left[\mbox{}_{2}^{}{{P}}{_{{k}-{m},{l}-{n}}^{mnkl}}
({p}^{2}{mk}+{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})+
\mbox{}_{2}^{}{{P}}{_{{k}+{m},{l}-{n}}^{mnkl}}({p}^{2}{mk}
-{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})\right.\\
\left.+\,\mbox{}_{2}^{}{{P}}{_{{k}-{m},{l}+{n}}^{mnkl}}
({p}^{2}{mk}-{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})+
\mbox{}_{2}^{}{{P}}{_{{k}+{m},{l}+{n}}^{mnkl}}
({p}^{2}{mk}+{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})\right]
\end{array}\!\right\}\nonumber\\
&&\,\,\,\,-\,\frac{\mu }{1-\mu }\frac{1}{64}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\vspace*{6pt}\nonumber\\
&&\,\,\,\,\times\, \frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
\left\{\!\!\begin{array}{@{}l@{}}
\left[\mbox{}_{2}^{}{{M}}{_{{k}-{m},{l}-{n}}^{mnkl}}({p}^{2}{mk}+
{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})+
\mbox{}_{2}^{}{{M}}{_{{k}+{m},{l}-{n}}^{mnkl}}
({p}^{2}{mk}-{q}^{2}{nl}+\alpha _{{mn}}\alpha _{{kl}})\right.\\
\left.+\,\mbox{}_{2}^{}{{M}}{_{{k}-{m},{l}+{n}}^{mnkl}}({p}^{2}{mk}
-{q}^{2}{nl}-\alpha _{{mn}}\alpha _{{kl}})+\mbox{}_{2}^{}
{{M}}{_{{k}+{m},{l}+{n}}^{mnkl}}({p}^{2}{mk}+{q}^{2}{nl}+
\alpha _{{mn}}\alpha _{{kl}})\right]
\end{array}\!\!\right\}\nonumber\\
\label{eqA6}
\end{eqnarray}}\unskip

$\mbox{}_{2}^{}{{P}}{_{rs}^{mnkl}}$ and
$\mbox{}_{2}^{}{{M}}{_{rs}^{mnkl}}$ being Fourier coefficients
defined by
{\begin{eqnarray}
{sh}(\alpha _{{mn}}+\alpha _{{kl}})
(\eta -{d}_{2})&=&\sum_{{u}=0}^{{\infty }}
\sum_{{v}=0}^{{\infty }}\Delta _{{u}0}
\Delta _{{v}0}\,\mbox{}_{2}^{}{{P}}{_{uv}^{mnkl}} 
\cos uX \cos vY\label{eqA7}\Seqnsplit
{sh}(\alpha _{{mn}}-\alpha _{{kl}})
(\eta -{d}_{2})&=&\sum_{{u}=0}^{{\infty }}
\sum_{{v}=0}^{{\infty }}\Delta _{{u}0}
\Delta _{{v}0}\,\mbox{}_{2}^{}{{M}}{_{uv}^{mnkl}} 
\cos uX \cos vY.
\label{eqA8}
\end{eqnarray}}\unskip

Using expressions~(\ref{eqA3}) and~(\ref{eqA6}), the kinetic energy 
can be expressed in a compact form by
{\begin{eqnarray}
\mathrm{KE}&=&\frac{1}{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times 
\frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
[\mbox{}_{1}^{}{\Omega }{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}\neq 
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
[\mbox{}_{1}^{}{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}^{2}(\alpha _{{mn}}{d}_{1})}\times 
[\mbox{}_{1}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
[\mbox{}_{2}^{}{\Omega }{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}
\neq \alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
[\mbox{}_{2}^{}{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}^{2}(\alpha _{{mn}}{d}_{2})}\times 
[\mbox{}_{2}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]
\label{eqA9}
\end{eqnarray}}\unskip
where we have
{\begin{eqnarray*}
\mbox{}[{\Theta }_{mnkl}^{+}]&=&\left[\begin{array}{l}
p^{2}mk+q^{2}nl+\alpha _{mn}\alpha _{kl}\\
p^{2}mk-q^{2}nl+\alpha _{mn}\alpha _{kl}\\
p^{2}mk-q^{2}nl-\alpha _{mn}\alpha _{kl}\\
p^{2}mk+q^{2}nl-\alpha _{mn}\alpha _{kl}
\end{array}\right]\quad  [{\Theta }_{mnkl}^{-}]=
\left[\begin{array}{l}
p^{2}mk+q^{2}nl-\alpha _{mn}\alpha _{kl}\\
p^{2}mk-q^{2}nl-\alpha _{mn}\alpha _{kl}\\
p^{2}mk-q^{2}nl+\alpha _{mn}\alpha _{kl}\\
p^{2}mk+q^{2}nl+\alpha _{mn}\alpha _{kl}
\end{array}\right]\\ 
\mbox{}[\mbox{}_{i}^{}{\Omega }{_{mnkl}^{+}}]&=&
\left[\begin{array}{l}
\mbox{}_{i}^{}{{P}}{_{{k}-{m},{l}-{n}}^{mnkl}}\vspace*{3pt}\\
\mbox{}_{i}^{}{{P}}{_{{k}-{m},{l}+{n}}^{mnkl}}\vspace*{3pt}\\
\mbox{}_{i}^{}{{P}}{_{{k}+{m},{l}-{n}}^{mnkl}}\vspace*{3pt}\\
\mbox{}_{i}^{}{{P}}{_{{k}+{m},{l}+{n}}^{mnkl}}
\end{array}\right]\quad (i=1,2) \quad [\mbox{}_{i}^{} {\Omega }{_{mnkl}^{-}}]
=\left[\begin{array}{l}
\mbox{}_{i}^{}{{M}}{_{{k}-{m},{l}-{n}}^{mnkl}}\vspace*{3pt}\\
\mbox{}_{i}^{}{{M}}{_{{k}-{m},{l}+{n}}^{mnkl}}\vspace*{3pt}\\
\mbox{}_{i}^{}{{M}}{_{{k}+{m},{l}-{n}}^{mnkl}}\vspace*{3pt}\\
\mbox{}_{i}^{}{{M}}{_{{k}+{m},{l}+{n}}^{mnkl}}
\end{array}\right]\quad (i=1,2)\\
\mbox{}[\mbox{}_{i}^{}{\Lambda }{_{mn}^{}}]&=&
\left[\begin{array}{l}
{a}_{{k}-{m},{l}-{n}}+4(-1)^{{i}}{d}_{{i}}\\
{a}_{{k}-{m},{l}+{n}}\\
{a}_{{k}+{m},{l}-{n}}\\
{a}_{{k}+{m},{l}+{n}}
\end{array}\right].
\end{eqnarray*}}\unskip

\subsection{The term of potential energy}\label{secA2}

In the right-hand side of~(\ref{eq7}), the averaged potential energy,
${V}_{{g}}$, due to gravity only is defined as 
{\begin{equation}\label{eqA10}
{V}_{{g}}=\frac{1}{4\uppi ^{2}}\int_{0}^{2\uppi }
\int_{0}^{2\uppi } \left[\frac{1}{2}\eta ^{2}\right]
\mathrm{d}X \mathrm{d}Y=\frac{1}{8}\sum_{{m}=0}^{\infty }
\sum_{{n}=0}^{\infty }\Delta _{{m}0}
\Delta _{{n}0}{{a}}_{{mn}}^{2}.
\end{equation}}\unskip

On the other hand, for the term of the averaged capillary energy
defined by 
{\begin{equation}\label{eqA11}
{V}_{\tau }=\frac{\kappa }{4\uppi ^{2}}\int_{0}^{2\uppi }
\int_{0}^{2\uppi }\left(\sqrt{1+{p}^{2}
\left(\frac{\partial \eta }{\partial {X}}\right)^{2}
+{q}^{2}\left(\frac{\partial \eta }{\partial {Y}}
\right)^{2}}-1\right)\mathrm{d}X\mathrm{d}Y
\end{equation}}\unskip
we introduce the following Fourier expansion 
{\begin{equation}\label{eqA12}
\sqrt{1+p^{2}\left(\frac{\partial \eta }{\partial X}\right)^{2}+
q^{2}\left(\frac{\partial \eta }{\partial Y}\right)^{2}}=
\sum_{m=0}^{\infty }\sum_{n=0}^{\infty }
\Delta _{m0}\Delta _{n0}R_{mn}\cos mX\cos nY.
\end{equation}}\unskip

That enables us to compute easily the integral to obtain then 
{\begin{equation}\label{eqA13}
{V}_{\tau }=\frac{\kappa }{4}({R}_{00}-4).
\end{equation}}\unskip

$R_{00}$ is computed by using the fast Fourier transform (FFT)
algorithm.         

\subsection{The term of averaged temporal rate of variation of the
velocity potentials}\label{secA3}

The remaining terms represents the total temporal rate of variation of
the potentials and has the following expression
{\begin{equation}\label{eqA14}
\overline{\Phi}_{{t}}=-\frac{1}{1-\mu }
\frac{1}{4\uppi ^{2}}{{\int }}_{0}^{2\uppi }{{\int }}_{0}^{2\uppi }
{{\int }}_{-{d}_{1}}^{\eta }\varpi \,\Phi _{1{X}}
\mathrm{d}Z \mathrm{d}X \mathrm{d}Y-
\frac{\mu }{1-\mu }\frac{1}{4\uppi ^{2}}{{\int }}_{0}^{2\uppi }
{{\int }}_{0}^{2\uppi }{{\int }}_{\eta }^{{d}_{2}}
\varpi \,\Phi _{1{X}}\mathrm{d}Z \mathrm{d}X \mathrm{d}Y.
\end{equation}}\unskip

After substituting the expansions~(\ref{eq12}) and~(\ref{eq13}) 
into~(\ref{eqA14}), and after integration with respect to $Z$, one
obtains
{\begin{eqnarray}
\overline{\Phi }_{{t}}&=&-\frac{1}{1-\mu }\frac{\varpi}
{4\uppi ^{2}}\sum_{{m}=1}^{\infty }\sum_{{n}=0}^{\infty }{m}
\frac{{b}_{{mn}}}{\alpha _{{mn}}}\int_{0}^{2\uppi }
\int_{0}^{2\uppi }\frac{{sh}(\alpha _{{mn}}
(\eta +{d}_{1}))}{{ch}
(\alpha _{{mn}}{d}_{1})}\cos mX \cos nY
\mathrm{d}X \mathrm{d}Y\nonumber\\
&&+\,\frac{\mu }{1-\mu }\frac{\varpi }{4\uppi ^{2}}
\sum_{{m}=1}^{\infty }\sum_{{n}=0}^{\infty }{m}
\frac{{c}_{{mn}}}{\alpha _{{mn}}}\int_{0}^{2\uppi }
\int_{0}^{2\uppi }\frac{{sh}(\alpha _{{mn}}
(\eta -{d}_{2}))}{{ch}
(\alpha _{{mn}}{d}_{2})}\cos mX \cos nY 
\mathrm{d}X \mathrm{d}Y.\label{eqA15}
\end{eqnarray} }\unskip

The integrals of the right-hand-side suggest to use the coefficients of
Fourier's expansion of ${sh}(\alpha _{{mn}}(\eta +{d}_{1}))$ and
${sh}(\alpha _{{mn}}(\eta -{d}_{2}))$, that is;
{\begin{eqnarray*}
{sh}(\alpha _{{mn}}(\eta +{d}_{1}))
&=&\sum_{{i}=0}^{{\infty }}\sum_{{j}=0}^{{\infty }}
\Delta _{{i}0}\Delta _{{j}0}
(\mbox{}_{1}^{}{\omega }{_{ij}^{mn}})
\cos iX \cos jY\\
{sh}(\alpha _{{mn}}(\eta -{d}_{2}))
&=&\sum_{{i}=0}^{{\infty }}\sum_{{j}=0}^{{\infty }}
\Delta _{{i}0}\Delta _{{j}0}(\mbox{}_{2}^{}{\omega }{_{ij}^{mn}})
\cos iX \cos jY.
\end{eqnarray*}}\unskip

In doing so, $\overline{\Phi }_{{t}}$ is easily found as
{\begin{equation}\label{eqA16}
\overline{\Phi }_{{t}}=-\frac{1}{1-\mu }\frac{\varpi }
{4}\sum_{{m}=1}^{{\infty }}\sum_{{n}=0}^{{\infty }}
{m}\frac{{b}_{{mn}}}{\alpha _{{mn}}{ch}(\alpha _{{mn}}
{d}_{1})}\mbox{}_{1}^{}{\omega }{_{mn}^{mn}}+
\frac{\mu }{1-\mu }\frac{\varpi }{4}\sum_{{m}=1}^{{\infty }}
\sum_{{n}=0}^{{\infty }}{m}\frac{{c}_{{mn}}}
{\alpha _{{mn}}{ch}(\alpha _{{mn}}{d}_{2})}
\mbox{}_{2}^{}{\omega }{_{mn}^{mn}}.
\end{equation}}\unskip

Using expressions~(\ref{eqA9}),~(\ref{eqA10}),~(\ref{eqA13}),
and~(\ref{eqA16}), the averaged Lagrangian function is obtained as a
function of the coefficients and can be written in the form 
{\begin{eqnarray}
\overline{{L}}&=&\frac{1}{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times 
\frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
[\mbox{}_{1}^{}{\Omega }{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}
\neq \alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{1}){ch}
(\alpha _{{kl}}{d}_{1})}\times \frac{1}
{(\alpha _{{mn}}-\alpha _{{kl}})}
[\mbox{}_{1}^{}{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&+\,\frac{1}{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{b}_{{mn}}{b}_{{kl}}}
{{ch}^{2}(\alpha _{{mn}}{d}_{1})}\times 
[\mbox{}_{1}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}+\alpha _{{kl}})}
[\mbox{}_{2}^{}{\Omega }{_{mnkl}^{+}}]^{ {t}} 
[{\Theta }_{{mnkl}}^{+}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}
\neq \alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}(\alpha _{{mn}}{d}_{2}){ch}
(\alpha _{{kl}}{d}_{2})}\times 
\frac{1}{(\alpha _{{mn}}-\alpha _{{kl}})}
[\mbox{}_{2}^{}{\Omega }{_{mnkl}^{-}}] ^{{t} }
[{\Theta }_{{mnkl}}^{-}]\nonumber\\
&&-\,\frac{\mu }{64(1-\mu)}\underset{\alpha _{{mn}}=
\alpha _{{kl}}}{\underbrace{\sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }\sum_{{k}=1}^{\infty }
\sum_{{l}=0}^{\infty }}}\frac{{c}_{{mn}}{c}_{{kl}}}
{{ch}^{2}(\alpha _{{mn}}{d}_{2})}\times 
[\mbox{}_{2}^{}{\Lambda }{_{mn}^{}}]^{ {t}} 
[{\Theta }_{{mnmn}}^{+}]\nonumber\\
&&+\,\frac{1}{8}\sum_{{m}=0}^{\infty }\sum_{{n}=0}^{\infty }
\Delta _{{m}0}\Delta _{{n}0}{{a}}_{{mn}}^{2}
-\frac{1}{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{m}\frac{{b}_{{mn}}}{\alpha _{{mn}}{ch}
(\alpha _{{mn}}{d}_{1})}\mbox{}_{1}^{}
{\omega }{_{mn}^{mn}}\nonumber\\
&&+\,\frac{\mu }{4(1-\mu)}\varpi \sum_{{m}=1}^{\infty }
\sum_{{n}=0}^{\infty }{m}\frac{{c}_{{mn}}}{\alpha _{{mn}}{ch}
(\alpha _{{mn}}{d}_{2})}\mbox{}_{2}^{}
{\omega }{_{mn}^{mn}}+\frac{\kappa }{4}({R}_{00}-4).
\label{eqA17}
\end{eqnarray}}\unskip

The details of first and second derivatives of the averaged
lagrangian are available from the authors. They are useful to
calculate the terms of the variational equations~(\ref{eq18})
and~(\ref{eq19}) and to solve them.

\bibliographystyle{crunsrt}
\bibliography{crmeca20230171}
\refinput{crmeca20230171-reference.tex} 

\end{document}
