Bifurcation Analysis and Propagation Conditions of Free-Surface Waves in Incompressible Viscous Fluids of Finite Depth

Abstract

Viscous linear surface waves are studied at arbitrary wavelength, layer thickness, viscosity, and surface tension. We find that in shallow enough fluids no surface waves can propagate. This layer thickness is determined for some fluids, water, glycerin, and mercury. Even in any thicker fluid layers, propagation of very short and very long waves is forbidden. When wave propagation is possible, only a single propagating mode exists for a given horizontal wave number. In contrast, there are two types of non-propagating modes. One kind of them exists at all wavelength and material parameters, and there are infinitely many such modes for a given wave number, distinguished by their decay rates. The other kind of non-propagating mode that is less attenuated may appear in zero, one, or two specimens. We notice the presence of two length scales as material parameters, one related to viscosity and the other to surface tension. We consider possible modes for a given material on the parameter plane layer thickness versus wave number and discuss bifurcations among different mode types. Motion of surface particles and time evolution of surface elevation is also studied at various parameters in glycerin, and a great variety of behaviour is found, including counterclockwise surface particle motion and negative group velocity in wave propagation.

1. Introduction

Although capillary gravity surface waves have been studied for a long time, they are still a fascinating and interesting subject [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. The dispersion of linear and nonlinear waves, wave damping in deep water owing to viscosity, and weakly nonlinear waves in shallow water are classic problems that have long been studied [16]. Comprehending the dissipation and dynamics of internal waves occurring at the boundary between viscous fluids holds significant importance. Surfactant, pollutant, and fluid film impacts at interfaces are now hot subjects [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35]. Current research on surface waves remains ongoing, with a predominant emphasis on multi-layer models [18], damping rate [2,5,7,17,18,19,20,36,37], and experimental verifications [17,20]. The majority of research endeavours in this particular domain are founded upon the utilization of viscoelastic fluids [17,21]. Hence, the constancy of surface tension is no longer maintained. The Marangoni number has been identified as a significant factor in the determination of a fluid’s surface behaviour, owing to its ability to characterize the fluid’s tendency to flow as a result of surface tension gradients [20,22,23]. The aforementioned phenomenon exhibits potential practical implications in the field of hydrodynamics, particularly in microfluidics and the study of nanoscopic capillary waves. This assertion is supported by various scholarly sources [24,25,26,27,28]. Whilst recent theoretical and experimental investigations [29,30] have primarily focused on fluid in deep layers, it is important to note that the structure and velocity of surface waves are significantly impacted by finite depth.

It is noteworthy that the uncomplicated scenario of linear surface waves in a viscous fluid possessing a clean surface and finite depth has not yet been fully explored. Only few works have been identified that address this problem, as documented in the references [33,36,38,39]. References [36,38] exclusively examined the scenario where surface tension was absent. The reference [39] incorporated surface tension, although the dispersion relation it obtained did not align with those presented in references [36,38] under the condition of zero surface tension. The appendix of reference [33] examines linear surface waves in finite depth, neglecting surface tension, and providing an approximation that is valid only under weakly damped conditions. Additionally, it should be noted that the dispersion relations derived in the works [36,38] exhibit dissimilarities.

The current investigation involves the derivation of the dispersion relation for surface waves that are both viscous and gravity-capillary in nature. This is achieved for fluid depths that are arbitrary. The outcome of our study in the limit of zero surface tension concurs with the dispersion relation derived in Hunt’s work [38]. As far as current understanding allows, the present work contains the first correct dispersion relation for viscous surface waves, which also incorporates the influence of surface tension and finite fluid depth. The present study demonstrates concurrence with prior research [17,21] in the context of the deep water regime. Upon examining the implications of the dispersion relation, it becomes apparent that a multitude of over-damped modes exist across all parameter values in an infinite fashion. There exist two distinct types of over-damped modes, characterized by either monotonic or oscillatory behaviour along the vertical axis. The latter category is a perpetual presence (in actuality, there are an infinite number of such modes at a specified wave number), whereas there may be zero, one, or two modes of the former category at a given wave number.

It has been observed that a single propagating mode can exist in a specific direction. This mode experiences a bifurcation at both long and short wavelengths, resulting in the emergence of two over-damped modes. The parameter space wave number versus layer thickness (we actually use the parameter p, which is related to layer thickness, cf. Equation (23)) can be thus partitioned according to the presence or absence of propagating modes and types and number of over-damped modes. The present study includes a discussion, primarily through numerical demonstration, of the bifurcations that occur at the boundaries of the distinct parameter space regions.

This paper is organized in the subsequent fashion: In Section 2, a derivation and discussion of the dispersion relation is presented, which serves as the basis for our subsequent investigations. Section 3 of the paper investigates various modes, their dependence on parameters, and bifurcations. Section 4 is dedicated to analysing the minimum layer thickness required for wave propagation. Section 5 delves into the examination of particle motion at the surface, and Section 6 scrutinizes the temporal progression of surface elevations. Section 7 of the paper offers a comprehensive analysis of the findings.

2. Dispersion Relation

The propagation of surface waves in deep water has been studied by several works, with results that generally agree with each other [17,21,24,40,41]. They are based on linearization of the Navier–Stokes equation. Here, we also consider linear surface waves on a viscous, incompressible fluid layer of finite, constant depth h. The coordinates x and y are horizontal, z is vertical. The origin lies at the undisturbed fluid surface. Suppose that the flow corresponding to the surface wave does not depend on y.

The linearized Navier–Stokes equation may be written as

$$\begin{array}{c}\hfill \frac{\partial \mathbf{V}}{\partial t}-\nu ▵\mathbf{V}=\nabla \left(-\frac{P}{\rho }-gz\right).\end{array}$$

(1)

where g indicates acceleration due to gravity, $\mathbf{V}$ is velocity vector, P is pressure, $\nu$ is viscosity, and $\rho$ is fluid density. Since the right hand side is a full gradient, we have

$$\begin{array}{c}\hfill \frac{\partial }{\partial z}\left(\frac{\partial u}{\partial t}-\nu ▵u\right)=\frac{\partial }{\partial x}\left(\frac{\partial w}{\partial t}-\nu ▵w\right).\end{array}$$

(2)

Due to linearity, we may search for the horizontal u and vertical w velocity components in the following form:

$$\begin{array}{c}\hfill u(x,z,t)=f^{\prime }\left(z\right)e^{i(kx-\omega t)}\end{array}$$

(3)

$$\begin{array}{c}\hfill w(x,z,t)=-ikf\left(z\right)e^{i(kx-\omega t)}\end{array}$$

(4)

here the prime denotes derivative with respect to the argument. Further, k is a real, positive wave number, while $\omega$ is in general complex frequency, its imaginary part describing the damping. The specified x and t dependence is due to translational symmetry in horizontal direction and time, respectively, while the specific form of z dependence is forced by the incompressibility condition $\nabla ·\mathbf{V}=0$. Note that the stream function is given by $fexp\left(i\right(kx-\omega t\left)\right)$.

Putting Ansatz (3), (4) into Equation (2), we obtain

$$\begin{array}{c}\hfill f^{″″}+\left(i\frac{\omega }{\nu }-2k^2\right)f^{″}-\left(i\frac{\omega }{\nu }-k^2\right)k^{2}f=0.\end{array}$$

(5)

Equation (5) has exponential solutions $f=exp\left(\kappa z\right)$. For the exponent $\kappa$ we get

$$\begin{array}{c}\hfill {\kappa }^4+\left(i\frac{\omega }{\nu }-2k^2\right){\kappa }^2-\left(i\frac{\omega }{\nu }-k^2\right)k^{2}f=0.\end{array}$$

(6)

The solutions are

$$\begin{array}{c}\hfill {\kappa }_{1,2}=\pm k,\end{array}$$

(7)

$$\begin{array}{c}\hfill {\kappa }_{3,4}=\pm \sqrt{k^2-i\frac{\omega }{\nu }}.\end{array}$$

(8)

For brevity, we shall use the notation $\kappa$ for ${\kappa }_3=-{\kappa }_4$ and k for ${\kappa }_1=-{\kappa }_2$. The general solution for f may be given as

$$\begin{array}{c}\hfill f=a_{1}cosh\left[k(z+h)\right]+a_{2}sinh\left[k(z+h)\right]+b_{1}cosh\left[\kappa (z+h)\right]+b_{2}sinh\left[\kappa (z+h)\right],\end{array}$$

(9)

where $a_1$, $a_2$, $b_1$, $b_2$ are the integration constants and h stands for the fluid depth. Then boundary conditions at the bottom,

$$\begin{array}{c}\hfill u(z=-h)=w(z=-h)=0,\end{array}$$

(10)

imply

$$\begin{array}{c}\hfill a_1=A,b_2=B,b_1=-A,a_2=-\frac{\kappa }{k}B,\end{array}$$

(11)

expressed in terms of the new constants A and B. Hence for f we get

$$\begin{array}{c}\hfill f=Acosh\left[k(z+h)\right]-\frac{\kappa }{k}Bsinh\left[k(z+h)\right]-Acosh\left[\kappa (z+h)\right]+Bsinh\left[\kappa (z+h)\right].\end{array}$$

(12)

Upon integrating the x component of the Navier–Stokes equation with respect to x, we get the pressure as

$$\begin{array}{c}\hfill P=P_o-\rho gz-\rho e^{i(kx-\omega t)}\left(-\frac{\omega }{k}{f}^{\prime }-i\nu k{f}^{\prime }+i\frac{\nu }{k}{f}^{‴}\right).\end{array}$$

(13)

At the fluid surface we have zero boundary conditions, therefore (in linear approximation) we have

$$\begin{array}{c}\hfill \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}=0\end{array}$$

(14)

for the shear and

$$\begin{array}{c}\hfill P-2\rho \nu \frac{\partial w}{\partial z}=P_0-\sigma \frac{{\partial }^2\eta }{\partial x^2}\end{array}$$

(15)

for the pressure. Here $\eta =\eta (x,t)$ stands for the deviation of the fluid surface from equilibrium and $\sigma$ is the surface tension.

Equation (14) implies

$$\begin{array}{c}\hfill f^{″}+k^{2}f=0\end{array}$$

(16)

at $z=0$, while Equation (15) implies

$$\begin{array}{c}\hfill -g\eta +\frac{\sigma }{\rho }\frac{{\partial }^2\eta }{\partial x^2}-e^{i(kx-\omega t)}\left(-\frac{\omega }{k}{f}^{\prime }-i\nu k{f}^{\prime }+i\frac{\nu }{k}{f}^{‴}\right)+2i\nu k{f}^{\prime }{e}^{i(kx-\omega t)}=0.\end{array}$$

(17)

In linear approximation, we have at the surface

$$\begin{array}{c}\hfill \frac{\partial \eta }{\partial t}=w.\end{array}$$

(18)

Note that on the right hand side we may set $z=0$. Putting here the expression of w (i.e., Equation (4)) and combining the result with Equation (17), in terms of $\kappa$ (cf. Equation (8)), we have

$$\begin{array}{c}\hfill \left(1+\frac{\sigma }{g\rho }{k}^2\right)k^{2}f+\frac{{\nu }^2}{g}({\kappa }^2-k^2)({\kappa }^2+2k^2)f^{\prime }-\frac{{\nu }^2}{g}({\kappa }^2-k^2)f^{‴}=0.\end{array}$$

(19)

Here, again, $z=0$. By now inserting the solution (12) into Equations (16) and (19) we obtain a linear homogeneous system of equation for the quantities A and B. A non-trivial solution can be found only if the determinant of the $2\times 2$ matrix formed out of the coefficients of A and B is zero. The resulting dispersion relation, expressed in terms of dimensionless quantities, sounds

$$\begin{array}{cc}& K\left(QsinhKcoshQ-KcoshKsinhQ\right)(1+s{K}^2)+p\left[-4K^{2}Q\left(K^2+Q^2\right)\right.\\ & \left.+Q\left(Q^4+2K^2{Q}^2+5K^4\right)coshKcoshQ-K\left(Q^4+6K^2{Q}^2+K^4\right)sinhKsinhQ\right]=0\end{array}$$

(20)

Here

$$\begin{array}{c}\hfill K=kh\end{array}$$

(21)

$$\begin{array}{c}\hfill Q=\kappa h\end{array}$$

(22)

$$\begin{array}{c}\hfill p=\frac{{\nu }^2}{g{h}^3}\end{array}$$

(23)

$$\begin{array}{c}\hfill s=\frac{\sigma }{\rho g{h}^2}\end{array}$$

(24)

Given parameters p and s, and scaled wave number K, a solution Q of Equation (20) yields the angular frequency (cf. Equations (8) and (22))

$$\begin{array}{c}\hfill \omega =-i\frac{\nu }{h^2}\left(K^2-Q^2\right).\end{array}$$

(25)

For the ratio of the coefficients A and B we get (cf. Equations (12) and (16))

$$\begin{array}{c}\hfill \frac{B}{A}=\frac{2K^{2}cosh\left(K\right)-(K^2+Q^2)cosh\left(Q\right)}{2QKsinh\left(K\right)-(K^2+Q^2)sinh\left(Q\right)}.\end{array}$$

(26)

Let us compare dispersion relation (20) with those available in the literature. The dispersion relation (20) for $\sigma =0$ coincides with the result of Hunt [38], albeit he used different notations. Meur [36] derived the dispersion relation

$$\begin{array}{ccc}& & 4K^{2}Q(K^2+Q^2)+4Q{K}^3[Qsinh\left(K\right)sinh\left(Q\right)-Ksinh\left(K\right)cosh\left(Q\right)]\hfill \\ & & -{(K^2+Q^2)}^2[Qcosh\left(K\right)cosh\left(Q\right)-Ksinh\left(K\right)sinh\left(Q\right)]\hfill \\ & & -KR{e}^2[Qsinh\left(K\right)sinh\left(Q\right)-Kcosh\left(K\right)sinh\left(Q\right)]=0\hfill \end{array}$$

(27)

for $\sigma =0$. Here $Re=\frac{hc}{\nu }$ where $c=\sqrt{gh}$, hence, according to Equation (23),

$$p=\frac{1}{R{e}^2}.$$

(28)

Note that $Re$ looks formally like a Reynolds number, but it has a quite different physical meaning. Reynolds number traditionally means the ratio of the nonlinear term to the viscous term in the Navier–Stokes equation. Now, there is no nonlinear term present. Further, in the definition of the present $Re$, the velocity c is not the fluid velocity but the propagation velocity of waves in shallow fluid.

Equation (27) is almost identical with (20), but there is a difference in the the bracket of $KR{e}^2[\dots ]$. The first term should be $Qsinh\left(K\right)cosh\left(Q\right)$. Sanochkin [39] also introduced a dispersion relation. All terms in [39] are the same as Equation (20) if $cosh\left(a\sqrt{q}\right)$ and $sinh\left(\frac{a}{\sqrt{1+q}}\right)$ change to $cosh\left(a\sqrt{1+q}\right)$ and $sinh\left(a\sqrt{1+q}\right)$, respectively. Note that q in his notation represents $\frac{Q^2}{K^2}-1$ and $a\equiv K$.

In the limit of $K\to \infty$, the dispersion relation (20) agrees with [21]. The author studied the dispersion relation in terms of dimensionless variables in three layers.

$$K\left(1+s{K}^2\right)+p\left(Q-K\right)\left(Q^3+Q^{2}K+3K^{2}Q-K^3\right)=0$$

(29)

Parameters p and s may be expressed in terms of the viscous length scale and the capillary length scale, which is related to surface tension, respectively. These parameters characterize the physical mechanisms that govern the oscillatory motion of capillary waves.

$$\begin{array}{c}\hfill {\mathit{\ell }}_{\nu }={\left(\frac{{\nu }^2}{g}\right)}^{1/3}\end{array}$$

(30)

$$\begin{array}{c}\hfill {\mathit{\ell }}_{\sigma }={\left(\frac{\sigma }{\rho g}\right)}^{1/2},\end{array}$$

(31)

namely,

$$\begin{array}{c}\hfill p={\left(\frac{{\mathit{\ell }}_{\nu }}{h}\right)}^3\end{array}$$

(32)

$$\begin{array}{c}\hfill s={\left(\frac{{\mathit{\ell }}_{\sigma }}{h}\right)}^2.\end{array}$$

(33)

It is obvious that the ratio

$$\begin{array}{c}\hfill \mu =\frac{{\mathit{\ell }}_{\sigma }}{{\mathit{\ell }}_{\nu }}=\frac{s^{1/2}}{p^{1/3}}\end{array}$$

(34)

is a material parameter (apart from g), and does not depend on the layer thickness. On the other hand, when parameter dependence is studied, often this is done by changing the layer thickness of a given material. In that case it is advisable to use parameters $\mu$ and p, while

$$\begin{array}{c}\hfill s={\mu }^2{p}^{2/3}.\end{array}$$

(35)

Rajan [21] also utilizes dimensionless variables in dispersion relation that involve the properties of the fluid and of the interface. We are interested in comparing wave-Reynolds ${\sigma }_{Rajan}$ number with s and p

$${\sigma }_{Rajan}=\frac{(1-R)+s{K}^2}{(1+R)p{K}^3}$$

(36)

where $R=\frac{{\rho }_{upper}}{{\rho }_{lower}}$. The density ratio R is restricted to take values less than or equal to unity due to the requirement of static stability. As $R\to 1$ in the fluid system, the effects of gravity tend to vanish [21]. In the case $R=0$, Equation (36) reduces to the well-known form of free fluid surface if the Marangoni number is zero. Therefore, we have

$${\sigma }_{Rajan}=\frac{1+s{K}^2}{p{K}^3}$$

(37)

In this study the Marangoni number P equals to zero, which represents zero interfacial elasticity and implies that the interface is clean.

Note that if Q is a solution of Equation (20), then so is $-Q$. On the other hand, this sign does not matter when calculating $\omega$ or f (cf. Equation (12)). Henceforth we assume that the real part of Q is positive, and thus $tanhQ\to 1$ when $\left|Q\right|\to \infty$.

The dispersion Equation (20) is valid for waves in a general system of fluids with arbitrary density, viscosity, surface tension, and depth. Here we modify Equation (20) for small viscosity and large wavelength cases that are applicable to specific systems. We also obtain numerical solutions. To this end, Equation (20) was expanded as a 50-order polynomial and all roots were found. The non-physical roots arising from the expansion were eliminated by the requirements that $\omega$ satisfy the original equation. Our analyses reveal the presence of bifurcation in some physical solutions, which are discussed in Section 2.4.

2.1. Small Viscosity Case

In the small viscosity case, Equation (20) may be solved approximately. In that case $p\to 0$ and $\left|Q\right|\to \infty$. This implies that in leading order Equation (20) reduces to

$$\begin{array}{c}\hfill KtanhK(1+s{K}^2)+p{Q}_0^4=0,\end{array}$$

(38)

or

$$\begin{array}{c}\hfill Q_0^2=-i\sqrt{\frac{KtanhK(1+s{K}^2)}{p}}.\end{array}$$

(39)

Here the negative sign has been chosen in order to get a positive real part of angular frequency via Equation (25). Further, according to the convention mentioned above, we have

$$\begin{array}{c}\hfill Q_0=\frac{1-i}{\sqrt{2}}{\left(\frac{KtanhK(1+s{K}^2)}{p}\right)}^{1/4}.\end{array}$$

(40)

A systematic expansion in terms of $p^{1/4}$ leads in the next two orders to [38]

$$\begin{array}{c}\hfill Q=Q_0-\frac{K}{2sinh\left(2K\right)}-\frac{K^2}{2Q_0}\frac{Y^2+6Y+5}{Y(Y+4)}.\end{array}$$

(41)

Here $Y=4{sinh}^{2}K$. To this order we have for the angular frequency

$$\begin{array}{ccc}\hfill \omega & =& \left[\sqrt{\left(gk+\frac{\sigma k^3}{\rho }\right)tanh\left(kh\right)}-\frac{\sqrt{2\nu k^2\sqrt{\left(gk+\frac{\sigma k^3}{\rho }\right)tanh\left(kh\right)}}}{2sinh\left(2kh\right)}\right]\hfill \\ & -i& \left[\frac{\sqrt{2\nu k^2\sqrt{\left(gk+\frac{\sigma k^3}{\rho }\right)tanh\left(kh\right)}}}{2sinh\left(2kh\right)}+2\nu k^2\frac{Y^2+5Y+2}{Y(Y+4)}\right].\hfill \end{array}$$

(42)

The first term of the real part is the well-known dispersion relation of surface waves in ideal fluids. As for damping, the leading term is the first one in the second bracket, proportional to $\sqrt{\nu }$, except in deep fluid. In deep fluid ($K\to \infty$) this term vanishes and one gets the well-known damping exponent $2\nu k^2$ [17,19]. The result (42) was first published (for $\sigma =0$) in reference [42] and then to higher orders in reference [38]. Note that taking into account surface tension is formally equivalent with replacing p with $p/(1+s{K}^2)$.

2.2. Behaviour at Large Wavelengths

At $K=0$ Equation (20) reduces to

$$\begin{array}{c}\hfill Q^{5}cosh\left(Q\right)=0\end{array}$$

(43)

(independently of p and s), which can be solved analytically: either

$$\begin{array}{c}\hfill Q=0\end{array}$$

(44)

or

$$\begin{array}{c}\hfill Q=i(2n+1)\frac{\pi }{2},n=0,1,2,\dots \end{array}$$

(45)

The corresponding frequencies are $\omega =0$ and

$$\begin{array}{c}\hfill \omega =-i\frac{\nu {\pi }^2}{4h^2}{(2n+1)}^2.\end{array}$$

(46)

Note that the order of the limits $\nu \to 0$ and $K\to 0$ does matter. If we take the limit $\nu \to 0$ first, we get ideal fluid, and if we take $K\to 0$ first, we get a limit where viscosity dominates and no wave propagation is possible.

2.3. Numerical Methods for Finding Roots

Equations (44) and (45) proved to be important technically, as solutions of Equation (20) could be obtained numerically from the differential equation (obtained directly by differentiating Equation (20) with respect to K)

$$\frac{dQ}{dK}=-\frac{\partial D/\partial K}{\partial D/\partial Q}$$

(47)

where D stands for the left hand side of Equation (20). Then Equations (44) and (45) play the role of initial conditions at $K=0$. At bifurcations we get singular behaviour, which is avoided on the complex K plane, afterwards we return to real K values. This method worked, and one could even choose different branches by going around the singularity from left or right on a half circle, yet it was somewhat in a state-of-art, as we had to experiment to get the correct radius of the half-circles. Another method we applied to get Q was direct numerical solution of Equation (20). In this case, bracketing solutions were a nontrivial task. Note that we compared the results obtained with these two different methods, and found an excellent agreement. A graphical solution was also helpful in locating roots.

2.4. Numerical Results for Frequencies

Having obtained Q for parameters p, s and K one may calculate mode frequencies. We present such results for a given material—glycerine—at some selected p values and as a function of scaled wave number K. The grand picture is this: at any given parameter settings one obtains infinitely many solutions, organized into branches as K changes (cf. Figure 1). In the following, the mode frequencies belonging to branch n will be denoted by ${\mathsf{\Omega }}^{\left(n\right)}\left(K\right)$. They emanate from values (46) or from zero. Most of them remain purely imaginary, but the lowest two may collide when increasing K. At such a collision, a bifurcation takes place, as in Figure 1 and Figure 2a, and two imaginary solutions Q may combine to complex solutions with non-zero real parts. Note that this bifurcation is also present in the parameters of Figure 3, but cannot be seen at the given resolution (it happens at very small wave numbers). The presence of complex Q solutions manifests itself in the appearance of a real part for the frequency. This means that wave propagation becomes possible. In that case the imaginary parts are the same, and the real parts differ in their signs only. As the value of p is increased, after a while such a collision no longer occurs, as in Figure 2c. This means that no wave propagation is possible because all modes are over-damped. The intermediate situation is approximately shown in Figure 2b. The mechanism of this possible bifurcation will be studied in more details in the next section (cf. Figure 4 and Figure 5). A bifurcation that is closely related to the one we find here at short wavelengths has already been found by Refs. [21,37]. In those works an infinitely deep fluid was studied.

In contrast, the bifurcation at long wavelength has been found explicitly here, and can be seen in Figure 1 and Figure 2a, where imaginary parts of the two lowest branches collide at a small wave number or long wavelength (we consider the wave number to be the control parameter when we speak about bifurcations). Below that wave number, all frequencies are purely imaginary, while above that wave number a real part arises (not shown in the figures), hence wave propagation becomes possible, as mentioned above. This bifurcation is also shown in Figure 6 on the $p-K$ parameter plane. There, wave propagation is possible in region IV, and a bifurcation at low wave numbers takes place when for a sufficiently low parameter p the red line between region I and region IV is crossed. For small wave numbers this line is parabolic, as quantified by Equation (75). This equation readily implies that the critical wavelength is ${\lambda }_c\propto h^{5/2}$, hence this bifurcation appears only in fluids of finite depth. The physical reason why this bifurcation happens can be understood from a simple order of magnitude estimate as follows. Let us consider a finite volume of the fluid, vertically from top to bottom, horizontally from a trough to the nearest crest. The main restoring force is gravity. If $\eta$ is the surface elevation, the pressure difference between the two vertical sides of the fluid volume is $\propto g\rho \eta$, and the restoring force (per unit perpendicular horizontal distance) is $\propto -g\rho h\eta$. Friction with the bottom is $\rho \nu du/dz\propto \rho \nu u/h$ per unit surface (the characteristic vertical length scale is h, cf. Equation (74)), or $\propto \lambda \rho \nu u/h$ for the whole segment (numerical factors are neglected), where $\lambda$ is the wavelength. The mass of the fluid segment is $\rho h\lambda$. Hence, Newton’s second law reads in this case as

$$\begin{array}{c}\hfill \rho h\lambda \dot{u}=-\lambda \rho \nu u/h-g\rho h\eta .\end{array}$$

(48)

Further, we have the kinematic condition at the surface, $w=\dot{\eta }$, and the condition of incompressibility, $\partial u/\partial x+\partial w/\partial z=0$. The partial derivatives may be estimated as $\partial u/\partial x\propto u/\lambda$ and $\partial w/\partial z=w/h$. By inserting these relations into Equation (48), we have

$$\begin{array}{c}\hfill \ddot{\eta }+\frac{\nu }{h^2}\dot{\eta }+\frac{gh}{{\lambda }^2}\eta =0.\end{array}$$

(49)

It is clear that for long waves friction prevails and the wave becomes over-damped. Indeed, the Ansatz $\eta ={\eta }_{0}exp(-i\omega t)$ leads to

$$\begin{array}{c}\hfill \omega =-i\frac{\nu }{2h^2}\pm \sqrt{\frac{gh}{{\lambda }^2}-\frac{{\nu }^2}{4h^4}}.\end{array}$$

(50)

This shows that for very long waves frequencies are purely imaginary. When wavelength $\lambda$ is decreased, at around

$$\begin{array}{c}\hfill {\lambda }_c\propto \frac{g^{1/2}{h}^{5/2}}{\nu }\end{array}$$

(51)

a non-zero real part of the frequencies emerges. Note that Equation (51), obtained from rough order of magnitude estimates for the bifurcation point, agrees well with Equation (75).

There are a number of works where the bifurcation phenomena in surface waves were studied [23,43,44,45,46,47]. All these works considered the infinite depth fluid case. They fall into three categories: (1) a non-viscous case where bifurcation is due to nonlinearity [45,46]. (2) The presence of a surfactant where an interplay between dilational and capillary waves takes place [21,23,47]. (3) The presence of elastic membrane at the surface of the ideal fluid [43,44].

A further work, not belonging to the above categories, contains an ingenious application of bifurcation theory in order to prove the existence of linear waves and find the dispersion relation [48]. The system considered is a finite depth ideal (non-viscous) fluid without surface tension. Laminar rotational flow is assumed with constant vorticity (i.e., horizontal velocity is changing linearly with depth), and linear surface gravity waves at a fixed wavelength are considered. Stationary solutions are sought for in a coordinate system moving horizontally with constant speed compared to the fluid. This speed is considered to be the control parameter (in fact, the equivalent mass flux is used). Now, if the speed of the coordinate system is not equal to the phase speed of the waves, only the trivial, zero amplitude stationary solution exists, since the linear wave is a travelling wave as seen from that coordinate system. If, however, the speed of the coordinate system is equal to the phase speed, a further, nontrivial stationary solution emerges, whose existence is proved with mathematical exactitude under quite general circumstances. This technique has also been applied to study linear capillary waves in the same system, by neglecting gravity [49].

3. Parameter Dependence

3.1. Wave Modes

Plotting the real and imaginary parts of Equation (20) on the complex Q plane one usually observes several intersections, i.e., roots. They can be real, purely imaginary, or complex with non-zero real and imaginary parts (henceforth we shall use the term “complex” as an equivalent with “complex with non-zero real part and non-zero imaginary part”). In the first two cases the angular frequency (25) is purely imaginary, so these modes decay exponentially with time. Propagating modes are only possible if Q is complex.

3.1.1. Modes with Real Q

In this case, Equation (25) implies that $Q<K$, since the decay rate cannot be negative. We have found numerically that at a given parameter settings (s, p, and K) there can exist zero, one, or two real modes (due to the symmetry of the solutions, we consider roots only in the first quadrant). There is always a trivial solution $Q=K$. Then, we get $f=0$ from Equation (12), so this solution is irrelevant.

In case of nontrivial real solutions, the decay rate is always smaller than $\nu k^2$. Velocity components may be calculated from Equations (3), (4), (12) and (26). Provided that coefficient A is real, u is real, too, while w is purely imaginary, so there is a 90° phase shift in their x-dependent oscillations.

3.1.2. Modes with Imaginary Q

In this case the decay rate is always larger than $\nu k^2$. Such modes exist at any parameter setting, moreover, there are infinitely many of them. Indeed, if Q is purely imaginary as $Q=i\beta$ and its modulus is large, Equation (20) reduces to

$$\begin{array}{c}\hfill p\left[Q^{5}coshKcoshQ-K{Q}^{4}sinhKsinhQ\right]=0,\end{array}$$

(52)

or, substituting $Q=i\beta$,

$$\begin{array}{c}\hfill \beta coshKcos\beta -KsinhKsin\beta =0.\end{array}$$

(53)

Now, if $\beta =2n\pi$ (n being an integer), the left hand side is positive, and if $\beta =2n\pi +\pi /2$, the left hand side is negative. Therefore, between these values there is a root for any (arbitrarily large) n. The distance between imaginary roots is approximately constant and independent of viscosity.

As before, the phase of velocity components does not change with depth, while there is a 90° phase shift between the x and z components. It is interesting that imaginary Q causes an oscillatory behaviour with depth.

3.1.3. Modes with Complex Q

On the basis of our numerical investigations, we believe that at a given parameter setting at most one such mode can exist. If there exists one, then at the same parameter setting no mode with real Q can exist. This time the phases of velocity components do change with depth. An oscillatory dependence on depth is in principle present, but much less pronounced than in the imaginary case.

3.2. Bifurcations

If, at the given s and p parameters one adjusts K, the type and number of solutions changes. This is shown in Figure 4a–i where $s=7.44$ and $p=2.0$. The left hand side of Equation (20) divided by $K^{2}Qcosh\left(K\right)cosh\left(Q\right)$ is plotted (this factor is chosen to avoid exponential growth.) on the complex Q plane, the zero level lines of the real part in blue, and those of the imaginary part in red. Hence, intersections of red and blue curves are solutions of Equation (20). As wave number K grows, these solutions are moving on that plain and display bifurcations. Since there are other (sometimes non-physical) solutions present, we denoted by black arrows the solutions we wish to focus on. In Figure 4a one can see a real solution. A bifurcation is shown in Figure 4b, where the real solution and its mirror image collide at the origin, and a further increase of the wave number gives rise to an imaginary solution (cf. Figure 4c, lower arrow) and its mirror image.

Increasing the scaled wave number K further, one can observe that that two imaginary solutions in Figure 4c collide (see Figure 4d) and give rise to a complex solution (see Figure 4e). This complex solution gradually goes down to the real axis (cf. Figure 4e–h) and decays to two real solutions (Figure 4i) which survive any further increase of K.

Choosing the parameter values $s=7.44$ and $p=3.0$, one has a different scenario, see Figure 5a–c. This time only a single bifurcation takes place, namely, an imaginary solution and its mirror image collide at the origin (see Figure 5b) and give rise to a (second) real solution (cf. Figure 5c).

These observations allow us to find the mathematical conditions of the bifurcations, which mean transitions among different types of possible modes. These conditions appear as curves on the K-p plane and partition that plane according to the types of possible solutions. Since the p-dependence is simpler than the K-dependence, we find these curves at constant K by varying p, rather than varying K at constant p, as has been done above when demonstrating the bifurcations. The curves are certainly the same.

Figure 4. Zero level lines of the real (blue) and imaginary (red) part of Equation (20) divided by $K^{2}Qcosh\left(K\right)cosh\left(Q\right)$ plotted on the complex Q plane at parameter values $p=2.0$ and $s=7.44$.

Figure 4. Zero level lines of the real (blue) and imaginary (red) part of Equation (20) divided by $K^{2}Qcosh\left(K\right)cosh\left(Q\right)$ plotted on the complex Q plane at parameter values $p=2.0$ and $s=7.44$.

Figure 5. Zero level lines of the real (blue) and imaginary (red) part of Equation (20) divided by $K^{2}Qcosh\left(K\right)cosh\left(Q\right)$ plotted on the complex Q plane at parameter values $p=3.0$ and $s=7.44$. (a) $K=2.900$. (b) $K=2.926$. (c) $K=3.00$.

Figure 5. Zero level lines of the real (blue) and imaginary (red) part of Equation (20) divided by $K^{2}Qcosh\left(K\right)cosh\left(Q\right)$ plotted on the complex Q plane at parameter values $p=3.0$ and $s=7.44$. (a) $K=2.900$. (b) $K=2.926$. (c) $K=3.00$.

In all the above numerically discussed cases (Figure 4 and Figure 5) at the birth of new solutions one can observe that by approaching the critical parameter the zero level lines of the real part of Equation (20) develop an edge, and at bifurcation they become (locally) two straight lines crossing each other on the real and/or the imaginary axis (Figure 4b,d,h and Figure 5b). This means that at the bifurcation point the derivative of the line is undetermined, i.e., it has a form $\frac{0}{0}$. Note in passing that it is equivalent with the condition that the second derivative becomes infinite. Let us denote for brevity the real part of Q with $\alpha$ and its imaginary part with $\beta$ (i.e., $Q=\alpha +i\beta$), further, the left hand side of Equation (20), divided by $K^{2}Qcosh\left(K\right)cosh\left(Q\right)$ be $D(\alpha +i\beta )$. Then the zero level line of the real part of D is given by

$$\begin{array}{c}\hfill \Re \left(D(\alpha +i\beta )\right)=0.\end{array}$$

(54)

Here ℜ stands for the real part. Taking its derivative with respect to $\alpha$, we get

$$\begin{array}{c}\hfill \Re \left(D^{\prime }(\alpha +i\beta )+i{D}^{\prime }(\alpha +i\beta )\frac{d\beta }{d\alpha }\right)=0,\end{array}$$

(55)

hence the derivative of the curve is written as

$$\begin{array}{c}\hfill \frac{d\beta }{d\alpha }=\frac{\Re \left(D^{\prime }(\alpha +i\beta )\right)}{\Im \left(D^{\prime }(\alpha +i\beta )\right)}.\end{array}$$

(56)

Here ℑ stands for the imaginary part (without i). It follows that at bifurcation

$$\begin{array}{c}\hfill D^{\prime }(\alpha +i\beta )=0\end{array}$$

(57)

must be satisfied, together with $D(\alpha +i\beta )=0$. Now it is easily seen that $D(\alpha +i\beta )$ is real along both the real and the imaginary axis, while $D^{\prime }(\alpha +i\beta )$ is real along the real axis and imaginary along the imaginary axis. Therefore, in case of a bifurcation on the real axis we have

$$\begin{array}{ccc}& & D\left(\alpha \right)=0\hfill \end{array}$$

(58)

$$\begin{array}{ccc}& & D^{\prime }\left(\alpha \right)=0,\hfill \end{array}$$

(59)

i.e., two real equations for the two real parameters $\alpha$ and p (at constant K and $\mu$). Similarly, in case of a bifurcation on the imaginary axis we have

$$\begin{array}{ccc}& & D\left(i\beta \right)=0\hfill \end{array}$$

(60)

$$\begin{array}{ccc}& & i{D}^{\prime }\left(i\beta \right)=0,\hfill \end{array}$$

(61)

again two real equations for the two real parameters $\beta$ and p.

Equation (20) shows that $D\left(Q\right)$ may be expressed as

$$\begin{array}{c}\hfill D\left(Q\right)=\left(1+s{K}^2\right)F\left(Q\right)+pG\left(Q\right)\end{array}$$

(62)

where

$$\begin{array}{cc}\hfill F\left(Q\right)& =\frac{tanhK}{K}-\frac{tanhQ}{Q},\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill G\left(Q\right)=& -4\frac{K^2+Q^2}{coshKcoshQ}+\left(\frac{Q^4}{K^2}+2Q^2+5K^2\right)\hfill \\ \hfill & -\left(Q^4+6K^2{Q}^2+K^4\right)\frac{tanhK}{K}\frac{tanhQ}{Q}.\hfill \end{array}$$

(64)

In terms of these functions we have (cf. Equations (58) and (59))

$$\begin{array}{cc}\hfill & F\left(\alpha \right)G^{\prime }\left(\alpha \right)-G\left(\alpha \right)F^{\prime }\left(\alpha \right)=0\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill & \frac{p}{1+s{K}^2}=-\frac{F\left(\alpha \right)}{G\left(\alpha \right)}\hfill \end{array}$$

(66)

and (cf. Equations (60) and (61))

$$\begin{array}{cc}\hfill & F\left(i\beta \right)G^{\prime }\left(i\beta \right)-G\left(i\beta \right)F^{\prime }\left(i\beta \right)=0\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill & \frac{p}{1+s{K}^2}=-\frac{F\left(i\beta \right)}{G\left(i\beta \right)}\hfill \end{array}$$

(68)

for bifurcations on the real and imaginary axis, respectively. In these cases, Equation (65) or Equation (67) is solved numerically to get $\alpha$ or $\beta$, respectively, then the solution is inserted into Equations (66) and (68), respectively, which in turn are solved for p at constant K and $\mu$, taking into account Equation (35).

As for the transformation of an imaginary solution to a real one at the origin, mentioned above, for the critical $p\left(K\right)$ line one obtains

$$\begin{array}{c}\hfill \frac{p}{1+s{K}^2}=-\frac{F\left(0\right)}{G\left(0\right)}=\frac{1-\frac{tanhK}{K}}{5K^2-4\frac{K^2}{cosh\left(K\right)}-K^{3}tanhK}.\end{array}$$

(69)

Since both $F\left(Q\right)$ and $G\left(Q\right)$ are even functions of Q, it follows that $F^{\prime }\left(0\right)=G^{\prime }\left(0\right)=0$, hence Equations (65) and (67) are automatically satisfied for $\alpha =0$ and $\beta =0$, respectively. The results are plotted in Figure 6a,b for glycerin and water, respectively. These two figures are very much different, not only in the shape of the curves, but in their vertical scales, in the first place. The values in Figure 6b are of order ${10}^9$. Such a huge difference needs an explanation. When we solve the Equations (66), (68), or (69) for p, we actually solve an equation

$$\begin{array}{c}\hfill \frac{p}{1+{\mu }^2{K}^2{p}^{2/3}}=a,\end{array}$$

(70)

where a is a function of K alone and is of order $0.1$. The solution of Equation (70) is

$$\begin{array}{c}\hfill p=\frac{1}{{\mu }^3{K}^3}{\mathsf{\Phi }}^{-1}\left({\mu }^3{K}^{3}a\right),\end{array}$$

(71)

where ${\mathsf{\Phi }}^{-1}\left(x\right)$ stands for the inverse of function $\mathsf{\Phi }\left(x\right)=x/(1+x^{2/3})$. It is easily seen that

$$\begin{array}{c}\hfill {\mathsf{\Phi }}^{-1}\left(x\right)\approx x,\mathrm{if}x\ll 1,\end{array}$$

(72)

and

$$\begin{array}{c}\hfill {\mathsf{\Phi }}^{-1}\left(x\right)\approx x^3,\mathrm{if}x\gg 1.\end{array}$$

(73)

Therefore, when ${\mu }^3{K}^{3}a\ll 1$, the solution is approximately $p=a$, as is the case for glycerin. In contrast, if ${\mu }^3{K}^{3}a\gg 1$, then the solution is $p={\mu }^6{K}^6{a}^3$. This happens in the case of water and mercury, and we get huge values for p, except for very small K values. These large p values are equivalent with extremely small layer widths (cf. Section 4). In fact, they are so small that the bifurcations in question are hardly observable in case of water or mercury as well as other fluids with low viscosity. They might be observable, however, in more viscous fluids such as glycerin.

In short, the physical reason of the difference between Figure 6a,b is that in case of water the viscosity is much smaller than in the case of glycerin (the surface tensions in the two cases are comparable) and this is reflected in the corresponding $\mu$ values ($58.40$ for water and $0.45$ for glycerin ). The two order of magnitude difference in the values of $\mu$ for water and glycerin is further magnified by the cubic functional dependence (73).

On the other hand, for small enough K values, condition ${\mu }^3{K}^{3}a\ll 1$ is satisfied, hence in that limiting case the parameter space regions become independent of $\mu$. This is natural, since in the long wavelength limit the impact of the surface tension vanishes.

In Figure 6 the red line means the onset of the creation of a complex solution at the imaginary line, i.e., the solution of Equations (67) and (68). The bue line is the same at the real line (cf. Equations (65) and (66)), while the black line is the onset of the crossover from the imaginary to the real solution at the origin (Equation (69)). It is obvious that these lines must have a common point. The lines partition the $K-p$ parameter space to four regions, denoted in the figures by Roman numbers:

I.

Only imaginary solutions (infinitely many of them) are present;

II.

There is a single real solution and there are infinitely many imaginary solutions;

III.

There are two real solutions and there are infinitely many imaginary solutions;

IV.

There is a single complex solution and there are infinitely many imaginary solutions.

Figure 6. The maximal parameters p versus the scaled wave-number K for glycerin (a) and water (b). Note the vertical scale in (b). The red line signifies the emergence of complex solutions at the imaginary line, while the blue line represents the real solutions obtained from Equations (65) and (66). The black line denotes the onset of the crossover from imaginary to real solutions at the origin.

Figure 6. The maximal parameters p versus the scaled wave-number K for glycerin (a) and water (b). Note the vertical scale in (b). The red line signifies the emergence of complex solutions at the imaginary line, while the blue line represents the real solutions obtained from Equations (65) and (66). The black line denotes the onset of the crossover from imaginary to real solutions at the origin.

Asymptotics and some special cases of the bifurcation curves are the following.

1.

Imaginary to complex Q (border between regions I and IV, the red curve in Figure 6):

For small K ($K\ll 1$) we have

$$\begin{array}{ccc}& & Q=(1.1127+0.2509K^2)i\hfill \end{array}$$

(74)

$$\begin{array}{ccc}& & p=0.53667K^2\left(1+(s-3.8674)K^2\right)\hfill \end{array}$$

(75)

2.

Complex to real Q (border between regions IV and III, the blue curve in Figure 6):

For large K we have

$$\begin{array}{ccc}& & Q=0.6823K\hfill \end{array}$$

(76)

$$\begin{array}{ccc}& & p=\frac{1.7200(1+s{K}^2)}{K^3}\approx \frac{1.7200s}{K}.\hfill \end{array}$$

(77)

This implies (cf. Equations (21)–(24) and Equations (30) and (31)) that the transition to the over-damped mode happens at wavelength

$$\begin{array}{c}\hfill \lambda =3.6530\frac{\rho {\nu }^2}{\sigma }=3.6530\frac{{\mathit{\ell }}_{\nu }^3}{{\mathit{\ell }}_{\sigma }^2}.\end{array}$$

(78)

It is remarkable that this value does not depend on layer thickness h or gravity of Earth g, as it is a purely material constant. Numerically, we get $\lambda =5\times {10}^{-8}\mathrm{m}$ for water, $7\times {10}^{-4}\mathrm{m}$ for glycerin and ${10}^{-9}\mathrm{m}$ for mercury. In our opinion, this effect cannot be observed in water or mercury, but might be observed in glycerin. It is questionable if hydrodynamics are applicable at the scales obtained for water and mercury, albeit at least one study states that it is applicable down to the nanoscales [28]. It is obvious that in case of these short waves gravity does not play a role. Here we demonstrate, that our more general setup leads to the right results in this limiting case. There are other situations, however, when gravity becomes important. This happens, e.g., at the bifurcation at long wavelengths. Gravity and surface tension become equally important when the two bifurcations (at large and small wavelength) are close to each other, as in Figure 2a.

3.

Imaginary to real Q (borders between regions I, II, and III, the black curve in Figure 6):

For small K we have

$$\begin{array}{c}\hfill p=\frac{1}{3}\left[1+\left(s-\frac{7}{5}\right)K^2\right].\end{array}$$

(79)

At $K\approx 4.9435$ parameter p diverges as

$$\begin{array}{c}\hfill p=\frac{0.8452+20.654s}{4.9435-K},\end{array}$$

(80)

For a given material this implies (cf. Equation (35))

$$\begin{array}{c}\hfill p=\frac{8.8\times {10}^3{\mu }^6}{{\left(4.9435-K\right)}^3}.\end{array}$$

(81)

4.

The common point of the bifurcation parameter curves (red, blue, and black lines in Figure 6 satisfies

$$\begin{array}{c}\hfill \frac{F\left(0\right)}{G\left(0\right)}=-\frac{p}{1+s{K}^2}=\underset{Q\to 0}{lim}\frac{F^{\prime }\left(Q\right)}{G^{\prime }\left(Q\right)}=\frac{F^{″}\left(0\right)}{G^{″}\left(0\right)},\end{array}$$

(82)

since $F^{\prime }\left(0\right)=G^{\prime }\left(0\right)=0$. This implies

$$\begin{array}{c}\hfill F\left(0\right)G^{\prime \prime }\left(0\right)-G\left(0\right)F^{\prime \prime }\left(0\right)=0.\end{array}$$

(83)

The solution of this equation yields for the coordinates of the common point

$$\begin{array}{ccc}& & K=2.4152\hfill \end{array}$$

(84)

$$\begin{array}{ccc}& & p=0.05307(1+5.8332s).\hfill \end{array}$$

(85)

4. Minimal Layer Thickness Necessary for Wave Propagation

Viscosity not only damps waves, but it can even prevent their propagation. Indeed, propagation, mathematically, a real part of the complex angular frequency, appears only in region IV (cf. Figure 6). This implies that no gravity-capillary waves can propagate if p is large enough (or, equivalently, if the layer width is small enough). Further, even if the layer thickness is larger than the critical value, neither very long, nor very short waves can propagate.

The critical layer thickness is found from the maximum point of the curves bordering region IV in Figure 6 (cf. Equations (65)–(68)). This depends on the material parameter $\mu$, so we present the results in

Table 1. Material parameters $l_{\nu }$, $l_{\sigma }$ and $\mu$ for some fluids.

Material	${\mathit{l}}_{\mathit{\nu }}$	${\mathit{l}}_{\mathit{\sigma }}$	$\mathit{\mu }$	${\mathit{p}}_{\mathbf{max}}$	h
water	$4.67\times {10}^{-5}$	$2.73\times {10}^{-3}$	58.40	$2.255\times {10}^9$	$3.561\times {10}^{-8}$
glycerin	$5.04\times {10}^{-3}$	$2.27\times {10}^{-3}$	0.45 1	0.085	$1.146\times {10}^{-2}$
mercury	$1.09\times {10}^{-5}$	$1.90\times {10}^{-3}$	174.86	$1.625\times {10}^{12}$	$9.271\times {10}^{-10}$

¹ Properties are at 20 °C.

. Clearly, for water and mercury the critical layer width is so extremely small that at such scales even the applicability of standard hydrodynamics is more than questionable.

5. Particle Motion at Surface

In this section, we focus on the numerical simulation of particle trajectories associated with wave patterns at the surface of the fluid. Since we are considering linear waves, both the horizontal and vertical motion is that of a damped harmonic oscillator. If we neglect damping, the resulting motion due to these two perpendicular harmonic oscillations results in a general elliptical motion. Note that in the presence of damping, the actual motion will be a spiral. The parameters of the ellipse (see Figure 7), main axis a, small axis b, and angle of the main axis made with x direction, $\phi$ are related to the horizontal and vertical amplitudes and phase difference. The relevant relations are the following:

$$\begin{array}{cc}\hfill & a^2=\frac{X_o^2{cos}^2\phi -Z_o^2{sin}^2\phi }{{cos}^2\phi -{sin}^2\phi },\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill & b^2=\frac{Z_o^2{cos}^2\phi -X_o^2{sin}^2\phi }{{cos}^2\phi -{sin}^2\phi }.\hfill \end{array}$$

(87)

Angle of main axis compared to horizontal, $\phi$, is given by

$$tan\phi =\frac{c^2-1+\sqrt{{(c^2+1)}^2-4c^2{sin}^2\gamma }}{2ccos\gamma }.$$

(88)

Here

$$c=\left|\frac{Z_o}{X_o}\right|$$

(89)

and

$$\gamma =\xi -\chi ,$$

(90)

where $\xi$ and $\chi$ are phases of the horizontal and the vertical oscillations, respectively:

$$\begin{array}{cc}\hfill & x=X_o{e}^{i\chi }·e^{-i\omega t}\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill & z=Z_o{e}^{i\xi }{e}^{-i\omega t}\hfill \end{array}$$

(92)

Explicitly,

$$\begin{array}{cc}\hfill X_o{e}^{i\chi }=& \frac{i}{\omega }K[Q(Q^2+3K^2)·(cosh\left(Q\right)cosh\left(K\right)-1)\hfill \\ \hfill & -K·(K^2+3Q^2)sinh\left(Q\right)sinh\left(K\right)]\hfill \end{array}$$

(93)

and

$$\begin{array}{cc}\hfill Z_o{e}^{i\xi }=& \frac{1}{\omega }K(Q^2-K^2)·[Kcosh\left(K\right)sinh\left(Q\right)\hfill \\ \hfill & -Q·cosh\left(Q\right)sinh\left(K\right)].\hfill \end{array}$$

(94)

According to Equations (93) and (94), when Q is a real, the value of $\gamma$ is $-\frac{\pi }{2}$, but if Q is purely imaginary, then the value of $\gamma$ is $\frac{\pi }{2}$.

These parameters are plotted versus scaled wave number K in Figure 8 at $p=0.001$. Wave propagation is permitted from $K=0$ to $K=13$, as seen in Figure 3 (Except for very small K values where the motion always becomes damped. At the resolution of the figure, this regime can not be seen). The motion is virtually linear ($\frac{b}{a}\ll 1$) for small values of K but quickly transforms into a circle with $\frac{b}{a}$ close to unity at $K\approx 1$. The angle that the main axis makes with the horizontal, $\phi$, is almost zero up to this point, then it rapidly takes 90 degree turn and remains around that. The ellipse gets narrower for wave numbers greater than unity ($K>1$) and at wave number $K=7.5$, it degenerates into a linear segment. The corresponding motion is nearly perpendicular to the surface. The ellipse then gets a little bit thicker. Both the horizontal and vertical motion are excessively over-dampened at $K\approx 13$. The change in the phase difference between horizontal and vertical motion is also displayed. It is clear that with K being approximately 1, the phase difference is approximately $\pi /2$. When the ellipse degenerates into a line segment at $K\approx 7.5$, the phase difference reaches $\pi$. At the $K\approx 13$ bifurcation point, the phase difference $\gamma$ hits $-\pi /2$.

6. Time Evolution of Surface Elevations

Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15 demonstrate the time evolution of a propagating wave in a fluid, along with the associated dispersion relations. We have chosen glycerin as the fluid medium due to its physical properties, as other fluids may require a thinner layer for observation of the phenomena. Then, the critical thickness for glycerin is 1.1 cm, for water is $3.6\times {10}^{-8}$ m, and for mercury is $9.3\times {10}^{-10}$ m, which is the reason we have chosen the glycerin for simulation. We limit our study to the lowest two branches of the dispersion relation, because, as shown in Figure 1 these branches have the longest lifetime and, i.e., they are the least damped.

We aim to establish a correlation between theoretical predictions and empirical observations. The findings indicate that wave propagation is detected at a particular wave-number of $K=1$, within a narrow range, when $p=0.077$, which is in proximity to the critical value of $p=0.086$. Conversely, no propagation is observed at longer wavelengths, such as $K=0.62$, or at shorter wavelengths, such as $K=2$, as illustrated in Figure 9, Figure 10 and Figure 11. Figure 9 depicts the over-damped mode characterized by a long wavelength. At a specific wavelength of $\lambda =12\mathrm{cm}$, a wavenumber of $K=0.62$, it can be observed that wave propagation is not feasible due to the nodes remaining in the unchanged position (cf. Figure 2a). It should be noted that surface tension plays a minor role in this context. The nodes in Figure 10 are observed to be in motion due to the occurrence of wave propagation at a wavenumber of $K=1$. This is in agreement with Figure 2a. At $K=2$, the wave in Figure 11 exhibits over-damping. Surface tension is a dominant phenomenon in this particular regime.

At $p=0.001$, the range of wavelengths that undergo propagation is wider. However, it is worth noting that a non-propagating long-wave mode develops at a wavelength of 8 m in a fluid layer of 5 cm, which poses a challenge for observation (refer to Figure 12). Figure 13 illustrates the observation of a travelling wave with a shorter wavelength of $\lambda =6.33\mathrm{m}$. The demonstration shows the limited wavenumber spectrum within which long wavelength over-damped modes exist. Figure 14 shows the propagation of waves at the same layer width, with a wavelength of $\lambda =3.16\mathrm{cm}$. It should be noted that the velocity is significantly greater, as evidenced by the time labels. The propagation of waves with a shorter wavelength, specifically those with a wavelength of $\lambda =2.11\mathrm{cm}$, is prohibited, as shown in Figure 15.

For the purpose of demonstrating the group velocities present in the propagating modes, we overlay two wavelengths that are in close proximity. Specifically, we consider the expression $sin(k_1(x-x_o)-{\mathsf{\Omega }}^{\left(1\right)}t)-sin(k_2(x-x_o)-{\mathsf{\Omega }}^{\left(2\right)}t)$, where both waves are propagating towards the right. Nevertheless, the direct observation of this phenomenon appears improbable due to significant damping. In order to show the motion of the envelope, the superposition is amplified at a rate of $e^{\left(\right|\gamma \left|t\right)}$, where $\gamma$ denotes the lesser decay rate. Figure 16 illustrates the superposition of two waves with two nearby wavenumbers ($K=0.75$ and $K=0.76$) at three consecutive time intervals. The direction of their travel is towards the right. At the onset, the envelope exhibits a node. At a subsequent time interval, the node expands to form a bottleneck. The position of the bottleneck is moving to the right, which indicates that the group velocity is also positive in this scenario. The phenomenon of broadening can be attributed to the disparate damping exhibited by two modes that are superimposed. Figure 17 describes the superimposition of two waves with closely spaced wavenumbers at the same fluid depth as previously. The present instance involves a reduction in wavelength, specifically with values of K equal to 0.94 and 0.95. Although the phase velocity retains a positive value, the position of the envelope remains unaltered in implying zero group velocity. The phenomenon of bottleneck widening persists due to the dissimilar damping exhibited by the two modes. Figure 18 illustrates the presence of two superimposed modes, which exhibit a comparatively reduced wavelength of $K=1.2$ and $K=1.21$. At present, although the phase velocity remains positive, the retrograde motion of the envelope indicates a negative group velocity. It should be noted that the frequency’s real part exhibits a negative slope versus K in this instance.

An arbitrary initial condition means specifying the velocity field at an instant of time everywhere within the fluid layer. In linear approximation one may decompose such an initial condition in terms of modes:

$$\begin{array}{c}\hfill \mathbf{V}(x,z,t=0)={\int }_{-\infty }^{\infty }dk\sum _{n=1}^{\infty }{a}_n\left(k\right){\mathbf{V}}_n(k,z)e^{ikx}\end{array}$$

(95)

Here $a_n\left(k\right)$ stands for suitable expansion coefficients, while

$$\begin{array}{c}\hfill {\mathbf{V}}_n(k,z)=\left(\begin{array}{c}Aksinh\left[k(z+h)\right]-B\kappa cosh\left[k(z+h)\right]-A\kappa sinh\left[\kappa (z+h)\right]+B\kappa cosh\left[\kappa (z+h)\right]\\ -iAkcosh\left[k(z+h)\right]+iB\kappa sinh\left[k(z+h)\right]+iAkcosh\left[\kappa (z+h)\right]-iBksinh\left[\kappa (z+h)\right]\end{array}\right)\end{array}$$

(96)

is the velocity field of a mode (cf. Equations (3), (4) and (12)). Note that $\kappa ={\kappa }_n\left(k\right)$ corresponds to the n-th solution of Equation (20), while the ratio of A and B is specified by Equation (26) (The quantities A and B are thus specified up to an arbitrary normalization factor), hence they also depend on $\kappa$, and thus on n and k.

If we take the Fourier transform of both sides of Equation (95) with respect to the x variable, we get

$$\begin{array}{c}\hfill \tilde{\mathbf{V}}(k,z,t=0)=\sum _{n=1}^{\infty }{a}_n\left(k\right){\mathbf{V}}_n(k,z),\end{array}$$

(97)

where

$$\begin{array}{c}\hfill \tilde{\mathbf{V}}(k,z,t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }dx{e}^{-ikx}\mathbf{V}(x,z,t).\end{array}$$

(98)

For arbitrary later times t we have

$$\begin{array}{c}\hfill \tilde{\mathbf{V}}(k,z,t)=\sum _{n=1}^{\infty }{a}_n\left(k\right){\mathbf{V}}_n(k,z)e^{-i{\mathsf{\Omega }}^{\left(n\right)}\left(k\right)t}.\end{array}$$

(99)

Most modes are strongly damped. Therefore, leaving them out of the decomposition (99) may not lead to a significant error, except initially for a very short time. If we keep only the lowest two branches, it is possible to formulate the initial value problem in terms of the surface profile and its time derivative. Note that in the range of wave numbers where propagation is possible, the two branches differ only in the sign of the real part of the frequency, allowing a description of both directions of propagation. Explicitly, we may formulate the initial value problem in wave number space (Fourier space) as follows. At a given wave number we have two modes, therefore the decomposition (99) reduces to

$$\tilde{\eta }(k,t)=c_1\left(k\right)e^{-i{\mathsf{\Omega }}^{\left(1\right)}t}+c_2\left(k\right)e^{-i{\mathsf{\Omega }}^{\left(2\right)}t}$$

(100)

for surface elevation at time t, where (cf. Equation (18))

$$c_n\left(k\right)=a_n\left(k\right)\frac{k}{{\mathsf{\Omega }}^{\left(n\right)}}\left\{Acosh\left[kh\right]-B\frac{\kappa }{k}sinh\left[kh\right]-Acosh\left[\kappa h\right]+Bsinh\left[\kappa h\right]\right\},$$

(101)

and

$$\begin{array}{c}\hfill \tilde{\eta }(k,t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }dx{e}^{-ikx}\eta (x,t).\end{array}$$

(102)

If the initial surface profile is $\tilde{\eta }(k,t=0)$, then

$$\tilde{\eta }(k,t=0)=c_1+c_2$$

(103)

should hold. Similarly, given the initial vertical velocity profile ${\tilde{\eta }}_t(k,t=0)$ we have

$${\tilde{\eta }}_t(k,t=0)=-i{\mathsf{\Omega }}^{\left(1\right)}{c}_1-i{\mathsf{\Omega }}^{\left(2\right)}{c}_2$$

(104)

From this we get for the coefficients $c_1$ an $c_2$

$$\begin{array}{cc}\hfill c_1& =\frac{i{\tilde{\eta }}_t-{\mathsf{\Omega }}^{\left(2\right)}\tilde{\eta }}{{\mathsf{\Omega }}^{\left(1\right)}-{\mathsf{\Omega }}^{\left(1\right)}}\hfill \end{array}$$

(105)

$$\begin{array}{cc}\hfill c_2& =\frac{i{\tilde{\eta }}_t-{\mathsf{\Omega }}^{\left(1\right)}\tilde{\eta }}{{\mathsf{\Omega }}^{\left(2\right)}-{\mathsf{\Omega }}^{\left(1\right)}}\hfill \end{array}$$

(106)

This allows one to solve the initial value problem within the limits of the approximation sketched above. On the other hand, such an approximation is completely equivalent with a second order differential equation for the surface elevation

$$\frac{{\partial }^2\tilde{\eta }}{\partial t^2}+i({\mathsf{\Omega }}^{\left(1\right)}\left(k\right)+{\mathsf{\Omega }}^{\left(2\right)}\left(k\right))\frac{\partial \tilde{\eta }}{\partial t}-{\mathsf{\Omega }}^{\left(1\right)}\left(k\right){\mathsf{\Omega }}^{\left(2\right)}\left(k\right)\tilde{\eta }=0$$

(107)

As an application of Equations (100)–(107), a narrow initial Gaussian profile (blue line) encompassing both propagating modes and over-damped modes is chosen (cf. Figure 19). Time evolution can be seen on the displayed profiles with different colours.

It is seen that the two peaks develop and radiate off symmetrically due to the propagating modes. Since there is no initial velocity, the peaks are in the same weight. Gradually, the Gaussian wave gives way to a much broader Gaussian shape profile, composed of non-propagating over-damped modes (yellow line). It is worth mentioning that in this example over-damped modes have a much smaller damping rate than the propagating ones. Therefore, the contributions of the propagating modes, the peaks, decay, and the longer living over-damped mode content of the initial profile becomes visible.

7. Conclusions

Linear viscous capillary-gravity waves were studied in a channel of constant depth, without restricting the parameters. The resulting dispersion relation is published for the first time. We explored all the modes numerically. Modes were labelled by horizontal wave number K, a continuous parameter, and vertical wave number Q, a discrete, complex quantity. We found that there were always infinitely many non-propagating (over-damped) modes. Propagation can only occur in the two modes with the smallest decay. In a sufficiently thin layer no propagation occurs at all. When increasing layer thickness, a bifurcation occurs, which shows up in the plot of imaginary parts of frequencies, such as a collision of the lowest lying branches of modes. After that collision, at increasing wave numbers, one can observe a merging and a subsequent split of these branches. The wave number range of the merged section frequencies have a non-zero real part (propagation). We stress that even at those depths where propagation becomes possible, propagation at very low and very high wave numbers is still prohibited. The over-damping at large wave numbers is already known [21,37]. The over-damping at small wave numbers has been discussed in Section 2.4. This is due to the fact that by increasing the wavelength, the restoring force (due to gravity) does not change, but the friction increases significantly. We also determined the minimal layer thickness necessary for wave propagation. Further, we studied surface motion. Assuming a monochromatic wave propagating to the positive x direction, we found that a surface particle in a viscous fluid could rotate both clockwise or counterclockwise, depending on the wave number. We also demonstrated the propagation or non-propagation of waves in a few cases. In order to illustrate that both positive, zero, or even negative group velocities can occur, the beat of two nearby wave numbers was displayed at a few consecutive time instants. Finally, with the assumption that the effect of fast decaying high lying branches was negligible, we kept the lowest two modes and formulated the solution of the initial value problem of surface motion in wave number space. As an application, the time evolution of a narrow initial Gaussian surface elevation with zero velocity was studied and a radiation of propagating modes in the form of two oppositely travelling bumps was observed. A slowly decaying wide Gaussian was left behind, consisting of large wavelength non-propagating modes.

It is remarkable how many interesting phenomena occur in this simple setup. It is an intriguing question what happens in the presence of an elastic surface at finite fluid depth, since two propagation modes, a capillary mode and a dilational mode, exist in deep water [23].

It is a natural question whether and how these results can have practical applications. As we already noted, in water most of the phenomena discussed are unobservable, because its viscosity is small and surface tension is large. In thick oils, however, these phenomena might be observed. For actual applications (industrial, or perhaps biological) one has to consider thin layers, but then other surface effects, such as pollution or surface elasticity, should be included, as well. It is tempting to think of the effect of the surfactant because they clearly decrease surface tension and with that $\mu$. However, with surfactants a wealth of new phenomena appear, e.g., the possibility of a dilational wave due to the presence of surface elasticity and interplay between this propagating mode with capillary waves (cf. [21,23,47]). Such phenomena are certainly not covered by the current paper.

Another potential application might be in geology (lava flows, glaciers) with especially large viscosity or planetology with different g parameters. In the first case, however, many other effects can be important, too, so the results of the present work could be used as the order of magnitude estimates at best. In case of other planets, a larger g decreases $\mu$, since it is proportional to $g^{-1/6}$, but at the same time parameter p is decreased, so an even thinner fluid layer (or larger viscosity) is necessary to prevent wave propagation, than in Earth. However, the minimal layer thickness was determined numerically, so it is not straightforward to predict the effect of changing g. It may happen that under suitable planetary circumstances the minimal layer thickness for some fluid that is present in abundance will be significantly larger than that for water on Earth.