On the non-uniqueness of the kernel of the Zakharov equation in intermediate and shallow water: the connection with the Davey–Stewartson equation

Abstract

The Zakharov equation describes the evolution of weakly nonlinear surface gravity waves for arbitrary spectral shape. For deep-water waves, results from the Zakharov equation are well established. However, for two-dimensional propagation, in intermediate and shallow water, there are problems related to the treatment of apparent singularities in the contribution of the wave-induced set-up to the evolution of the surface gravity waves. More specifically, the kernel in the integral term is characterized by a regular and an apparent singular contribution. Here, we show that the Davey–Stewartson equation can be directly derived from the Zakharov equation, also in the shallow water limit. This result provides guidance on how to treat the singular contribution to the evolution of the action variable. A relevant result that is obtained is that the growth rate obtained from the stability analysis of a plane wave in shallow water does not depend on the singular part of the kernel of the Zakharov equation.

1. Introduction

Surface gravity waves are usually described in the context of the potential flow of an ideal fluid. As discovered by Zakharov (1968), the resulting nonlinear evolution equations are obtained from a Hamiltonian, which is the total energy $E$ of the fluid, while the appropriate canonical variables are the surface elevation $\eta (\boldsymbol {x},t)$ and the value $\psi$ of the potential $\phi$ at the surface, $\psi (\boldsymbol {x},t)= \phi (\boldsymbol {x},z=\eta,t))$.

In deep water the approach of Zakharov is well established. However, in shallow water there is a problem because it is not immediately evident how one should treat the apparent singularities related to the wave-induced set-up and current. For one-dimensional propagation, a solution to this problem was suggested by Janssen & Onorato (2007) by taking the limit of an infinitely long wave train for a single wave and agreement for surface elevation, wave-induced current and the nonlinear dispersion relation was found with the results of Whitham (1974). For the case of two-dimensional propagation, there is still discussion on how to proceed (see e.g. Stiassnie & Gramstad 2009).

To avoid ambiguity in the integrals, in Gramstad (2014), it was strongly emphasized to explicitly separate the wave-induced mean flow and set-up from the surface gravity oscillations, resulting in a rather involved version of the Zakharov equation. Here, we follow a slightly different route and show that such an approach is superfluous. Being aware of the work in Onorato et al. (2009), where it is shown that the Boussinesq equation, after removing non-resonant interactions, corresponds to the Zakharov equation in the shallow-water limit, one may expect that the Davey–Stewartson (DS) equation, describing a wave group, follows from the Zakharov equation in the narrow-band approximation. This is explicitly shown in this paper by extending the analysis in Janssen & Onorato (2007) to two horizontal dimensions. This is in keeping with results found by Stiassnie & Shemer (1984), who were probably the first to make this connection. Indeed, from the DS equation, an explicit expression for the nonlinear transfer function is obtained which compares favourably with the narrow-band limit of the nonlinear transfer according to the Zakharov equation. A number of properties of the resulting nonlinear transfer can then be studied. For example, as pointed out already by Herterich & Hasselmann (1980), the narrow-band limit of the nonlinear transfer is not unique as it depends on the order in which the limit of the vanishing of the modulation wavenumbers is taken. Nevertheless, we show that the non-uniqueness is resolved when one has knowledge of the initial condition of the two-dimensional part of the mean wave-induced current. In addition, we also perform a stability analysis of a plane wave and it is shown that, despite the non-uniqueness of the nonlinear transfer coefficient, the growth rate of the Benjamin–Feir instability in intermediate and shallow water is unique. Just recently, Pezzutto & Shrira (2023) discussed the ambiguity problem and, based on a limiting procedure, they claimed that, despite the fact that the coefficients are singular, the integrals are not. They used this result for computing the Stokes frequency correction and the growth rate in the modulational instability. Here, our prospective is different, and we do not make any attempt to compute directly the integrals. We show that, for a Stokes wave train, the growth rates are independent of such ambiguity, as this is connected to the initial condition on the mean flow; moreover, we are able to associate the singular part of the integral with the auxiliary variable $Q$ (related to the mean flow) in the DS equation. In agreement with Davey & Stewartson (1974), we find that the frequency correction in such limit can take an arbitrary constant value.

The content of this paper is as follows. In § 2, we start from the DS equation and we transform this coupled set of equations for the wave action variable and mean potential to Fourier space. Then, the Fourier transform of the mean potential may be eliminated in favour of the Fourier transform of the action variable, resulting in a standard four-wave interaction equation for the action variable. According to the DS equation, the nonlinear interaction coefficient is now explicitly known and the contribution of the wave-induced mean potential is found to depend on the modulation wavenumber, in such a way that, in the limit of vanishing modulation wavenumber, the limiting value of the transfer function for a homogeneous background solution is an arbitrary constant. In § 3, we start from the Zakharov equation and we determine the narrow-band approximation of the nonlinear interaction coefficient. The resulting approximation is in perfect agreement with the results from the Fourier transform of the DS equation. As a consequence, it is found that, in the narrow-band limit, both of the nonlinear transfers from the DS equation and the Zakharov equation are not unique. Nevertheless, in § 4, a unique answer for the growth rate of the Benjamin–Feir instability is found.

2. The Davey–Stewartson equation in Fourier space

The starting point is the DS equation, of which there are several equivalent forms available. Since the Zakharov equation is formulated in terms of the action variable, we will take the equation that is based on this variable. Hence, we introduce the action variable $a(\boldsymbol {k},t)$ defined in such a way that the surface elevation $\eta (\boldsymbol x,t)$ reads

(2.1)\begin{equation} \eta(\boldsymbol x,t)=\int {\rm d}\boldsymbol{k}\left(\frac{\omega_k}{2g}\right)^{1/2} [a(\boldsymbol{k},t)+a^*(-\boldsymbol{k},t)]\exp\left({\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}\right), \end{equation}

where $\omega _k=\omega (k)=\sqrt {g k \tanh (k h)}$, with $g$ acceleration of gravity, $k=|\boldsymbol {k}|=\sqrt {k_x^2+k_y^2}$ the magnitude of the wavenumber and $h$ the constant depth. Assuming a narrow-band process with carrier wavenumber $\boldsymbol {k}=(k_0,0)$ propagating in the $x$-direction, and defining the Fourier transform according to

(2.2)\begin{equation} A(\boldsymbol{x},t) = \int {\rm d}\boldsymbol{k} a(\boldsymbol{k},t)\exp\{{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}\}, \end{equation}

at leading order, one finds from (2.1)

(2.3)\begin{equation} \eta(\boldsymbol x,t)=\left(\frac{\omega_0}{2g}\right)^{1/2} [A(\boldsymbol{x},t) \exp({{\rm i}k_0x})+A^*(\boldsymbol{x},t)\exp({-{\rm i}k_0x})], \end{equation}

where $\omega _0=\omega (k_0)$. The DS equation for the complex envelope $A=A(\boldsymbol {x},t)$ and for the auxiliary variable $Q=Q(\boldsymbol {x},t)$, related to the mean motion of the fluid, then becomes

(2.4)\begin{equation} \left.\begin{gathered} {\rm i}\frac{\partial A}{\partial t}+\lambda \frac{\partial^2A}{\partial x^2}+\mu \frac{\partial^2A}{\partial y^2}=\nu|A|^2A+\nu_1QA\\ \lambda_1\frac{\partial^2 Q}{\partial x^2}+\mu_1 \frac{\partial^2 Q}{\partial y^2}=\nu_2\frac{\partial^2|A|^2}{\partial y^2} \end{gathered}\right\}, \end{equation}

where

(2.5)\begin{equation} \left.\begin{gathered} \lambda=\frac{1}{2}\left.\frac{\partial^2\omega}{\partial k_x^2}\right|_{(k_0,0)},\quad \mu=\frac{v_g}{2k_0},\\ \nu=k_0^3\left[\frac{9T_0^4-10T_0^2+9}{8T_0^3}-\frac{1}{k_0h} \left\{\frac{(2v_g-c_0/2)^2}{c_S^2-v_g^2}+1\right\}\right],\\ \nu_1=\frac{k_0^4}{2\omega_0v_g}\left\{2c_0+v_g(1-T_0^2)\right\},\quad \lambda_1=c_S^2-v_g^2,\quad \mu_1=c_S^2,\\ \nu_2=\frac{1}{2}\frac{c_0}{T_0}c_S^2v_g\frac{2c_0+v_g(1-T_0^2)}{c_S^2-v_g^2}. \end{gathered}\right\} \end{equation}

Furthermore, $T_0=\tanh (k_0 h)$, $c_0=\omega _0/k_0$ is the phase speed, $v_g=\partial \omega /\partial k_x|_{(k_0,0)}$ is the group speed and $c_S=\sqrt {gh}$ is the shallow-water speed. We remark that the one-dimensional form of the DS equation with $Q=0$ agrees with earlier work by Hasimoto & Ono (1972) and Whitham (1974). The two-dimensional version of the DS equation was already found and studied by Benney & Roskes (1969).

Using the Fourier transform with wave vector $\boldsymbol {\kappa }=\{l,m\}$, e.g.

(2.6)\begin{equation} q_1= \frac{1}{4{\rm \pi}^2}\int {\rm d}\kern0.7pt \boldsymbol{x} Q(\boldsymbol{x},t) \exp\{-{\rm i}\boldsymbol{\kappa}_1\boldsymbol{\cdot}\boldsymbol{x}\}, \end{equation}

we rewrite the DS equation as follows:

(2.7)\begin{equation} \left.\begin{gathered} {\rm i}\frac{\partial a_1}{\partial t}-(\lambda l_{1}^2 +\mu m_{1}^2)a_1=\nu \int{\rm d} {\boldsymbol \kappa}_{2,3,4}a_2^*a_3a_4\delta_{1+2-3-4} +\nu_1\int {\rm d}{\boldsymbol \kappa}_{3,i}q_i a_3\delta_{1-3-i} \\ (\lambda_1 l_i^2 +\mu_1 m_i^2) q_i= \nu_2\int {\rm d}{\boldsymbol \kappa}_{2,4}(m_4-m_2)^2 a_2^*a_4\delta_{i+2-4}. \end{gathered}\right\} \end{equation}

From the second equation in (2.7), one may express $q_i$ in terms of the action variable, i.e.

(2.8)\begin{equation} q_i=\nu_2\int {\rm d}{\boldsymbol \kappa}_{2,4}\frac{(m_4-m_2)^2}{\lambda_{1} l_i^2 +\mu_1 m_{i}^2} a_2^*a_4\delta_{i+2-4}. \end{equation}

Note that, in the limit of $l_i$ and $m_i$ going to zero, to get (2.8), we are dividing by zero; we believe that this is the origin of the singularities in the subsequent analysis. Substitution (2.8) in the first equation, gives the standard form of a four-wave interaction equation

(2.9)\begin{equation} {\rm i}\frac{\partial a_1}{\partial t}=\varOmega_1 a_1 +\int{\rm d}{\boldsymbol \kappa}_{2,3,4}G_{1,2,3,4}a_2^*a_3a_4 \delta_{1+2-3-4}, \end{equation}

with the dispersion relation

(2.10)\begin{equation} \varOmega_1 = \lambda l_{1}^2 +\mu m_{1}^2, \end{equation}

while the nonlinear transfer coefficient becomes

(2.11)\begin{equation} G_{1,2,3,4} = \nu+\nu_1\frac{\nu_2(m_{4}-m_{2})^2} {\lambda_1(l_{4}-l_{2})^2+\mu_1(m_{4}-m_{2})^2}. \end{equation}

For the equation to be Hamiltonian, the coefficient must satisfy a number of symmetries: using the condition $\boldsymbol {\kappa }_1+\boldsymbol {\kappa }_2=\boldsymbol {\kappa }_3+\boldsymbol {\kappa }_4$, one may also write $G_{1,2,3,4}$ in terms of the difference modulation vector $\boldsymbol {\kappa }_1-\boldsymbol {\kappa }_3$. Moreover, $G_{1,2,3,4}$ must be invariant with respect to the exchange of labels ‘3’ and ‘4’. Combining these properties and defining

(2.12)\begin{equation} E_{i,j} = \frac{1}{4}\frac{(m_i-m_j)^2}{\lambda_1(l_i-l_j)^2+\mu_1(m_i-m_j)^2}, \end{equation}

one obtains for $G$ the symmetric form

(2.13)\begin{equation} G_{1,2,3,4} = \nu+\nu_1[\nu_2(E_{1,3}+E_{2,3}+E_{4,2}+E_{4,1})] .\end{equation}

From (2.13), it is immediately evident that the effect of the mean flow on the evolution of the modulations is significant. In the function $E$ both denominator and numerator depend on the square of the difference between the modulation wavenumbers in such a way that, if the size of the modulation wavenumbers decreases, the impact of the function $E$, and hence the mean flow, does not diminish. This is in sharp contrast with the effect of mean flow on modulations in in infinite water depth as given by the two-dimensional Dysthe equation; indeed, in Dungey & Hui (1979), it has been shown show, and further discussed by Janssen (1983), that the deep-water wave-induced current is a higher-order effect.

It should be clear, see figure 1, that, at the origin the function, $E$ is not well defined. This is a potential problem because, when determining the growth rate of the Benjamin–Feir instability, we need to evaluate the nonlinear transfer function for vanishing modulation wavenumber difference, i.e. $l_i=l_j$ and $m_i=m_j$. To see this better, we consider $E$ as a function of $K_x=l_i-l_j$ and $K_y=m_i-m_j$, and set all the coefficients to 1, to get $E = {K_y^2}/(K_x^2+K_y^2)$. Now, in the limit of vanishing modulation wavenumber difference, i.e. $K_x \rightarrow 0, K_y\rightarrow 0$, one would expect to find a unique answer. However, this is not the case. Considering the straight line through the origin, $K_y=\alpha K_x$, one finds that $E=\alpha ^2/(1+\alpha ^2)$, so that the value of $E$ at the origin depends on the limiting procedure, i.e. under what angle the origin is approached. Despite this non-uniqueness, we will find that the evolution of the modulations, e.g. the linear growth rate of the modulational instability, is found to be unique.

Figure 1. Wave-induced current function $E= {K_y^2}/(K_x^2+K_y^2)$ where $K_x=l_i-l_j$ and $K_y=m_i-m_j$ as function of $K_x$ and $K_y$.

In contrast, for the deep-water Dysthe equation the effect of the wave-induced current has, after simplification, the form $E = {K_x^2}/(\sqrt {K_x^2+K_y^2})$ and this has at the origin the unique value of zero because the numerator vanishes more rapidly towards zero than the denominator. This is immediately seen by considering once more the straight line $K_x=\alpha K_y$ so that $E=\alpha ^2 K_x/\sqrt {1+\alpha ^2}$, which vanishes unambiguously for vanishing $K_x$.

3. The Zakharov equation and its connection with the DS equation

The Zakharov equation reads

(3.1)\begin{equation} {\rm i}\frac{\partial a_1}{\partial t}=\omega_1 a_1 + \int T_{1,2,3,4}a_2^*a_3a_4 \delta_{1+2-3-4}\,{\rm d}{\boldsymbol k}_{2,3,4}, \end{equation}

where the nonlinear interaction coefficient $T_{1,2,3,4}$ is given in the Appendix. We are interested in an envelope equation so that it is useful to rewrite the absolute value of the wavevector $k$ in the dispersion relation as

(3.2)\begin{equation} k = k_x\sqrt{1+\frac{k_y^2}{k_x^2}}. \end{equation}

Because in the DS equation $k_y^2/k_x^2\ll 1$, i.e. the wavenumber in the $y$-direction is much smaller than the one in the $x$-direction, we assume that there is a dominant wave propagating in the $x$ direction with wave vector components $(k_0,0)$, and we consider the narrow-band approximation

(3.3) \begin{equation} \left.\begin{gathered} k_x=k_0+ l\\ k_y= m. \end{gathered}\right\} \end{equation}

Expanding the dispersion relation for small $l$ and $m$, we obtain

(3.4)\begin{equation} \omega_k=\omega_0+ v_g l +(\lambda l^2+\mu m^2). \end{equation}

In the Zakharov equation the first two terms can be removed by a rotation, and we are left with the dispersion relation of the DS equation, see (2.10). We now need to calculate the narrow-band version of the coefficient $T_{1,2,3,4}$. Its explicit form is given in the Appendix, and it consists of two contributions. The first one is denoted by $T^{(R)}_{1,2,3,4}$ and represents the regular contribution. It can be verified (see e.g. Janssen & Onorato 2007) that the narrow-band version of the regular contribution becomes

(3.5)\begin{equation} T^{(R)}_{1,1,1,1}= k_0^3\frac{9T_0^4-10T_0^2+9}{8 T_0^3}. \end{equation}

However, there is also a contribution from the non-resonant triads. This contribution is denoted by $T^{(S)}_{1,2,3,4}$; these terms involve the product of the second-order coefficients $V^{(-)}$. The difficulty of the calculation comes from the terms that, in the limit of a long wave group, give rise to an apparent singularity. For example, the term that contains at the denominator $\omega _3+\omega _{1-3}-\omega _1$ may give rise to a singular behaviour when $\boldsymbol {k}_1\simeq \boldsymbol {k}_3$.

The calculation that follows is similar to the one performed by Janssen & Onorato (2007), but now extended to the case of two-dimensional propagation. We shall only give results explicitly for the product term of the triads $(1,3,1-3)$ and $(4,2,4-2)$ and the triads $(3,1,3-1)$ and $(2,4,2-4)$. The other contributions follow by interchanging the indices 1 and 2. We need to calculate terms like the following one (see (A3)):

(3.6)\begin{equation} {-}V^{(-)}_{1,3,1-3}V^{(-)}_{4,2,4-2}\left[\frac{1}{\omega_3+\omega_{1-3}-\omega_1}+ \frac{1}{\omega_2+\omega_{4-2}-\omega_4}\right]. \end{equation}

After some algebra one finds that (see also (23) in Janssen & Onorato 2007)

(3.7)\begin{equation} V^{(-)}_{1,3,1-3} \approx \frac{1}{4\sqrt{2}} \left(\frac{k_0^2-q_0^2} {\omega_0}\sqrt{g \omega_{1-3}}+2 k_0(l_1-l_3) \left(\frac{g}{\omega_{1-3}}\right)^{1/2}\right);\end{equation}

moreover, we need to compute terms such as

(3.8)\begin{equation} \frac{1}{\omega_3-\omega_1+\omega_{1-3}}. \end{equation}

To this end, it is straightforward to show that, to the leading non-trivial order, we have

(3.9) \begin{equation} \omega_{1-3}= c_s \varDelta_{1-3}, \end{equation}

with

(3.10)\begin{equation} \varDelta_{1-3}= \sqrt{(l_1-l_3)^2+(m_1-m_3)^2}, \end{equation}

while

(3.11)\begin{equation} \omega_1-\omega_3= v_g (l_1-l_3). \end{equation}

Furthermore, the wavenumber condition $\boldsymbol {k}_1+\boldsymbol {k}_2=\boldsymbol {k}_3+\boldsymbol {k}_4$ implies certain relations for the difference vectors, e.g. $l_1-l_3=l_4-l_2$ and $m_1-m_3=m_4-m_2$. As a consequence, one finds that $\omega _{1-3}=\omega _{4-2}$ while $\omega _1-\omega _3=\omega _4-\omega _2$. By interchanging 1 and 2, it is found that $\omega _{2-3}=\omega _{4-1}$ and that $\omega _2-\omega _3=\omega _4-\omega _1$. Similar relations exist for the parameter $\varDelta$, e.g. $\varDelta _{1-3}=\varDelta _{4-2}$ and $\varDelta _{2-3}=\varDelta _{4-1}$.

Putting together the first and the fourth term of $T^{(S)}_{1,2,3,4}$ denoted by $2F_{1,3}$, and the second and the third term given by $2F_{2,3}$, one then finds

(3.12)\begin{equation} T^{(S)}_{1,2,3,4} = 2F_{1,3}+2F_{2,3}, \end{equation}

with

(3.13)\begin{align} F_{i,j}=-\frac{1}{16}\frac{k_0^3 c_s^2}{c_s^2 \varDelta_{i-j}^2-v_g^2 (l_i-l_j)^2} \left[\frac{(1-T_0^2)^2 \varDelta_{i-j}^2}{T_0}+ \left(\frac{4(1-T_0^2)g v_g}{c_s^2 \omega_0}+\frac{4}{k_0 h}\right)(l_i-l_j)^2\right]. \end{align}

Furthermore, exploiting the above mentioned relations between the difference vectors, one obtains the following symmetrical form:

(3.14)\begin{equation} T^{(S)}_{1,2,3,4} = F_{1,3}+F_{2,3}+F_{4,2}+F_{4,1}, \end{equation}

so that the narrow-band version of the nonlinear interaction coefficient $T_{1,2,3,4}$ becomes

(3.15)\begin{equation} T_{1,2,3,4} \approx T^{(R)}_{1,1,1,1}+F_{1,3}+F_{2,3}+F_{4,2}+F_{4,1}, \end{equation}

where $T^{(R)}_{1,1,1,1}$ is given by (3.5). To show that the narrow-band version of the Zakharov interaction coefficient is identical to the corresponding one from the DS equation (3.3), we have to further work on the form of the function $F_{i,j}$. Straightforward algebraic manipulations show that this function is of the form

(3.16)\begin{equation} F_{i,j}= \alpha_1+\alpha_2 -4\mu_1\alpha_2E_{i,j},\quad {\rm with}\ E_{i,j}=\frac{1}{4} \frac{(m_i-m_j)^2}{\lambda_1(l_i-l_j)^2+ \mu_1(m_i-m_j)^2}, \end{equation}

where $E_{i,j}$ is identical to the form given in (2.12). Furthermore, the constants are

(3.17a,b)\begin{equation} \alpha_1=-\frac{1}{16}k_0^3 \frac{(1-T_0^2)^2}{T_0}\quad {\rm and}\quad \alpha_2=-\frac{1}{16}\frac{k_0^3}{T_0(c_s^2-v_g^2)}(2 c_0+v_g(1-T_0^2))^2. \end{equation}

Using the expressions for $\mu _1$ and $\alpha _2$, it is then straightforward to establish that the modulation dependent part of (2.13) matches the one in (3.16) since $-4\mu _1\alpha _2=\nu _1\nu _2$.

Now, it remains to establish whether the constant parts of (2.13) and (3.17a,b) agree. After some algebra, using (32)–(33) of Janssen & Onorato (2007), i.e.

(3.18)\begin{align} -\frac{1}{4}k_0^3\frac{c_s^2}{c_s^2-v_g^2}\left[\frac{(1-T_0^2)^2}{T_0}+ \frac{4(1-T_0^2)gv_g}{c_s^2\omega_0}+\frac{4 }{k_oh}\right] =-\frac{k_0^3}{k_0h}\left[\frac{(2v_g-c_0/2)^2}{c_s^2-v_g^2}+1\right], \end{align}

an alternative expression for $\alpha _2$ may be obtained which establishes a connection between $\alpha _2$ and $\alpha _1$, i.e.

(3.19)\begin{equation} \alpha_2=-\alpha_1 -\frac{1}{4}\frac{k_0^3}{k_0h}\left[\frac{(2v_g-c_0/2)^2}{c_s^2-v_g^2}+1\right]. \end{equation}

As a consequence, one may write for the narrow-band version of the nonlinear interaction coefficient

(3.20)\begin{equation} T_{1,2,3,4} \approx T_{1D}+\nu_1[\nu_2(E_{1,3}+E_{2,3}+E_{4,2}+E_{4,1})] ,\end{equation}

with

(3.21)\begin{align} T_{1D} = T^{(R)}_{1,1,1,1}+4(\alpha_1+\alpha_2)=k_0^3 \frac{9T_0^4-10T_0^2+9}{8 T_0^3}-\frac{k_0^3}{k_0h} \left[\frac{(2v_g-c_0/2)^2}{c_s^2-v_g^2}+1\right]. \end{align}

Hence, the term $T_{1D}$ agrees with the term $\nu$ in the expression for $G_{1,2,3,4}$ in (2.13), which gives the nonlinear transfer for the DS equation. Note that $T_{1D}$ is in agreement with Whitham (1974), and it gives the nonlinear correction to the dispersion relation for a uniform weakly nonlinear wave train and the correction due to the one-dimensional wave-induced current. The sum of these two terms vanishes for $k_0h=1.363$, which signals for one-dimensional modulations the transition from unstable to stable. However, as already shown by Hayes (1973) for two-dimensional modulations, there is still modulational instability for $k_0h<1.363$.

4. Linear stability analysis of a uniform wave train

The main topic of this section is to study the linear stability results of a uniform wave train, realizing that the nonlinear transfer coefficient is not unique in the limit of zero modulation wavenumber (see § 2). Because the growth rate will be seen to depend explicitly on the value of the $T_{1,2,3,4}$ at the carrier wavenumber $\boldsymbol {k}_0=(k_0,0)$ and at the sideband wavenumber $\boldsymbol {k}_0+\boldsymbol {\kappa }$, the main question to be answered is as follows: How does the mentioned non-uniqueness affect the stability results?

Let us now revisit the stability analysis of a plane wave. The nonlinear dispersion relation for surface gravity waves is obtained immediately from the Zakharov equation (3.1). Consider the case of a single wave

(4.1)\begin{equation} a(\boldsymbol{k},t)=\hat{a}(t)\delta(\boldsymbol{k}-\boldsymbol{k}_0), \end{equation}

then substitution of (4.1) into (3.1) gives

(4.2)\begin{equation} \frac{{\rm d} \hat{a}(t)}{{\rm d} t}=-{\rm i}T_{NL}|\hat{a}(t)|^2\hat{a}(t), \end{equation}

where $T_{NL}=T_{0,0,0,0}$. Equation (4.2) may be solved at once by writing

(4.3)\begin{equation} \hat{a}(t)=a_0 \exp({-{\rm i}\varOmega t}), \end{equation}

where $\varOmega$ denotes the correction of the dispersion relation due to nonlinearity. It is given by

(4.4)\begin{equation} \varOmega= T_{NL}|a_0|^2. \end{equation}

The dependence of the dispersion relation on the wave amplitudes is known to have a profound impact on the time evolution of a weakly nonlinear wave train, but here we just confine ourselves to discussing the short time behaviour of a nonlinear wave train by means of a linear stability analysis.

Following Crawford et al. (1981) (see also Yuen & Lake 1982; Krasitskii & Kalmykov 1993), we now address the question as to whether a weakly nonlinear wave train is stable or not in finite water depth (see also Brinch-Nielsen & Jonsson 1986; Gramstad & Trulsen 2011). To test the stability of a uniform wave train, we perturb it by a pair of sidebands with wavenumber $\boldsymbol {k}_{\pm }=\boldsymbol {k}_0 \pm \boldsymbol {\kappa }$ and amplitude $A_{\pm }(t)$, e.g.

(4.5)\begin{equation} a= A_0\delta(\boldsymbol{k}-\boldsymbol{k}_0)+A_+\delta(\boldsymbol{k}-\boldsymbol{k}_+)+A_-\delta(\boldsymbol{k}-\boldsymbol{k}_-). \end{equation}

Assuming that the sideband amplitudes are small compared with the amplitude $A_0$ of the carrier wave and neglecting the square of small quantities, the following evolution equations for $A_{\pm }$ are found from the Zakharov equation (3.1):

(4.6)\begin{equation} {\rm i}\frac{{\rm d}}{{\rm d} t}A_{{\pm}}=T_{{\pm},\mp}a_0^2A_{{\mp}}^*\exp[-{\rm i}(\Delta\omega+ 2T_{NL}a_0^2)t]+2T_{{\pm},\pm}a_0^2A_{{\pm}}, \end{equation}

where

(4.7)\begin{equation} \left.\begin{gathered} T_{{\pm},\pm}=T(\boldsymbol{k}_0 \pm \boldsymbol{\kappa},\boldsymbol{k}_0,\boldsymbol{k}_0,\boldsymbol{k}_0 \pm \boldsymbol{\kappa}),\\ T_{{\pm},\mp}=T(\boldsymbol{k}_0 \pm \boldsymbol{\kappa},\boldsymbol{k}_0 \mp \boldsymbol{\kappa},\boldsymbol{k}_0,\boldsymbol{k}_0),\\ T_{NL} = T(\boldsymbol{k}_0,\boldsymbol{k}_0,\boldsymbol{k}_0,\boldsymbol{k}_0), \\ \Delta\omega=2\omega(\boldsymbol{k}_0)-\omega(\boldsymbol{k}_0+\boldsymbol{\kappa})-\omega(\boldsymbol{k}_0-\boldsymbol{\kappa}), \end{gathered}\right\} \end{equation}

and $a_0$ is the same quantity as given in (4.4).

By means of the substitution

(4.8)\begin{equation} \left.\begin{gathered} A_+=\hat{A}_+\exp\left[-{\rm i}\left(\tfrac{1}{2}\Delta\omega+T_{NL}a_0^2\right)t-{\rm i}\varOmega t\right], \\ A_-^*=\hat{A}_-^*\exp\left[+{\rm i}\left(\tfrac{1}{2}\Delta\omega+T_{NL}a_0^2\right)t-{\rm i}\varOmega t\right], \end{gathered}\right\} \end{equation}

where $\varOmega$ is still unknown, a set of differential equations is obtained that contains no explicit time dependence. A non-trivial solution is then found provided $\varOmega$ satisfies the dispersion relation

(4.9) \begin{align} \varOmega &= (T_{+,+}-T_{-,-})a_0^2 \nonumber\\ &\quad \pm \bigg\{\!-T_{+,-}T_{-,+}a_0^4+\left[-\tfrac{1}{2}\Delta \omega+a_0^2(T_{+,+}+T_{-,-}-T_{NL})\right]^2\!\bigg\}^{{1}/{2}}. \end{align}

We have instability provided that the term under the square root is negative (see Crawford et al. 1981).

Let us now determine the dispersion relation explicitly for the narrow-band version of the nonlinear transfer coefficient given in (3.20). Evaluating the parameters given in (4.7) with $\boldsymbol {\kappa }=(l,m)$, one finds

(4.10)\begin{equation} \left.\begin{gathered} T_{+,+} =T_{-,-}=T_{1D}+\frac{1}{2} \nu_1\nu_2\frac{ m^2}{\lambda_1 l^2+\mu_1 m^2}+2\nu_1\nu_2E_0, \\ T_{+,-}=T_{-,+}=T_{1D}+\nu_1\nu_2\frac{ m^2}{\lambda_1 l^2+\mu_1 m^2}, \\ T_{NL}=T_{1D}+4\nu_1\nu_2E_0,\\ \Delta\omega=-2(\lambda l^2+\mu m^2). \end{gathered}\right\} \end{equation}

Here, $T_{1D}$ is the part of the nonlinear transfer coefficient that is independent of the modulation wavenumber $\boldsymbol {\kappa }=(l,m)$ and $E_0$ is the value of $E$ at the origin. Now, since $T_{+,+}=T_{-,-}$ and $T_{+,-}=T_{-,+}$, the dispersion relation (4.9) simplifies to

(4.11) \begin{equation} \varOmega={\pm} \big\{-T^2_{+,-}a_0^4 +\big[-\tfrac{1}{2}\Delta\omega+a_0^2(2T_{+,+}-T_{NL})\big]^2\big\}^{{1}/{2}}. \end{equation}

The most remarkable point to make now is that the term $(2T_{+,+}-T_{NL})$, upon using (4.10), does not depend on the undefined value $E_0$ as

(4.12)\begin{equation} 2T_{+,+}-T_{NL}=T_{1D}+\nu_1\nu_2\frac{m^2}{\lambda_1 l^2+\mu_1 m^2}=T_{+,-}.\end{equation}

Since also the term $T_{+,-}$ does not depend on the undefined value $E_0$, this implies that the dispersion relation for $\varOmega$ is well defined. Further simplification, using $T_{+,-}=2T_{+,+}-T_{NL}$, gives as dispersion relation

(4.13)\begin{equation} \varOmega^2=(\lambda l^2+\mu m^2)[2T_{+,-}a_0^2+(\lambda l^2+\mu m^2)], \end{equation}

which is in perfect agreement with stability results obtained directly from the DS equation. Hayes (1973) and Davey & Stewartson (1974) showed that the wave train is unstable if

(4.14)\begin{equation} T_{+,-}(\lambda l^2+\mu m^2) < 0. \end{equation}

Note that $\lambda = \omega _0''/2$ is always negative while $\mu =v_g/2k_0$ is positive, and $T_{1D}$ changes from negative to positive as $k_0h$ decreases beyond $k_0h=1.363$ (Hasimoto & Ono 1972). The two-dimensional effect involving the modulation wavenumbers $l$ and $m$ tends to enhance the Benjamin–Feir instability. Hayes (1973) has shown that it is always possible to choose $l$ and $m$ in such a way that criterion (4.14) is satisfied, but instability is practically non-existent for shallow-water waves in the range $0\leq k_0h\leq 0.5$.

Finally, it is remarked that, while $T_{+,-}$ is not unique at the origin, $\varOmega$ from (4.13) is well defined because $T_{+,-}$ is multiplied by $\Delta \omega$, which vanishes at the origin.

It is clear that the parameter $T_{+,-}$ plays an important role in the modulational instability. Therefore, in order to obtain confidence in our results, we compare this parameter with alternative approaches. From (4.10), we know that, in the narrow-band approximation, we have

(4.15)\begin{equation} T_{+,-}=T_{1D}+\nu_1\nu_2\frac{m^2}{\lambda_1 l^2+\mu_1 m^2}. \end{equation}

Some effort has been made to try to validate the result given in (4.15) using direct numerical computation of the original interaction coefficient. In figure 2, we show the ratio of the narrow-band approximation in (4.15) and the numerical result obtained by computing the element matrix of the Krasitskii's interaction coefficient. This last coefficient is evaluated following Janssen & Onorato (2007), where the wavenumbers have been perturbed by a small amount in such a way that the perturbations (with size of order $10^{-5}$) satisfy the wavenumber resonance conditions. Therefore, the computation of the last coefficient is also an approximation, and one might expect deviations much larger than machine precision. The analytical result was obtained by (3.20) where from the outset we used the wavenumber condition $\boldsymbol {k}_1+\boldsymbol {k}_2=\boldsymbol {k}_3+\boldsymbol {k}_4$ to simplify (3.20), so that it only depends on the first three wavenumbers, thus

(4.16)\begin{equation} T_{1,2,3} = T_{1D} + 2\nu_1\nu_2(E_{1,3} + E_{2,3}). \end{equation}

Figure 2. Dependence of $T_{1,2,3,4}$ on the dimensionless depth parameter $k_0h$. The case $k_1=k_0(\cos \theta,+\sin \theta )$, $k_2=k_0(\cos \theta,-\sin \theta )$, $k_3=k_4=k_0(\cos \theta,0)$ is chosen. Shown is the analytical result in (4.15) normalized with the numerical result for different values of $\theta$.

The conclusion, looking at figure 2, is that there is a reasonable agreement between the analytical result and the numerical implementation of the Krasitskii interaction coefficients. This implies that for numerical calculations of the nonlinear transfer in shallow and intermediate water, i.e. the procedure suggested by Janssen & Onorato (2007), seems to give valid results.

5. Discussion and conclusions

We have shown that in the narrow-band approximation the DS equation follows from the Zakharov equation. This means that the Zakharov equation contains all the physics presented by the DS equation. Moreover, we have shown that the growth rate of the two-dimensional version of the Benjamin–Feir instability obtained from the Zakharov equation agrees with the result from the DS equation. While it is shown that the growth rate does not depend on the undefined value $E_0$, this is an entirely different matter for the Stokes frequency correction given by (4.4). It is seen that this correction does depend explicitly on $T_{0,0,0,0}$, therefore there is an ambiguity, which can only be resolved by making a choice for $E_0$. To see this, we look for an homogeneous solution for the envelope $A({\boldsymbol x},t)$ of the DS equation in (2.4) by making the ansatz $A({\boldsymbol x},t)=A_0 \exp ({-\textrm {i}\varOmega t})$, with $A_0$ a constant; the second of the equations in (2.4) becomes

(5.1)\begin{equation} \lambda_1\frac{\partial^2 Q}{\partial x^2}+\mu_1\frac{\partial^2 Q}{\partial y^2}=0; \end{equation}

therefore, for a constant (in space) envelope, a solution of the above equation is that $Q$ takes any arbitrary constant, $Q= \textrm {const}=Q_0$. This is all in agreement with Davey & Stewartson (1974).

Hence, the nonlinear frequency shift depends on the chosen value of $Q_0$. We note that, if $Q$ is a constant then, using (2.6), $q_{\boldsymbol {\kappa }}=Q_0\delta (\boldsymbol {\kappa })$; therefore, with $A_0=\textrm {const}$, we get $A_{\boldsymbol {\kappa }}=A_0 \delta ({\boldsymbol \kappa })\exp ({-\textrm {i}\varOmega t})$, and the apparent singular part of the kernel in (3.20) of the Zakharov equation must be taken as a constant, i.e.

(5.2)\begin{equation} \lim_{\boldsymbol {\kappa}_1,\boldsymbol {\kappa}_2\boldsymbol {\kappa}_3,\boldsymbol {\kappa}_4\rightarrow 0} (E_{1,3}+E_{2,3}+E_{4,2}+E_{4,1})=Q_1, \end{equation}

with $Q_1=Q_0/(\nu _2|A_0|^2)$. The value of $Q_1$ can be also chosen as 0, so that the kernel becomes identical to the one-dimensional case computed in Janssen & Onorato (2007). We remark, once more, that the choice of the constant $Q_0$ is immaterial for the stability analysis of a plane wave.

Finally, the non-uniqueness of the interaction kernel gives indeed rise to uncertainties in simulations by means of, e.g. wave forecasting models. If one would just calculate the nonlinear transfer without any extra precautions, the calculation would explode at locations in the direction–wavenumber grid where the numerator and the denominator of the interaction kernel vanish. In order to avoid such explosions, a procedure was suggested in Janssen & Onorato (2007) in such a way that finite answers were obtained in agreement with known results from Whitham (1974) in one dimension, and in the present paper for two dimensions with the DS equation. The procedure is to perturb the relevant wavenumbers by a slight amount in such a way that the wavenumber resonance condition is satisfied. We underline that, in Yang, Yao & Zhang (2022), the growth rates of the modulational instability using the Janssen & Onorato (2007) approach, extended to two-dimensional propagation, have been compared with numerical results obtained from a higher order spectral method simulation of the Euler equations and good agreement of the maximum growth rates was reported.