**1. Introduction**

Waves on the surface of rivers, seas and oceans became one of the main objects of theoretical research as the first fundamental equations of continuum mechanics and closed systems of equations were formulated [1–3]. Due to a large difference in the physical properties of the atmosphere and water [4], the theory of waves evolved with the approximation of the homogeneity and immutability of the density of the contacting media. B. Franklin in the 18th century observed the water–olive oil interface motion in a ship's lighting lamp and pointed out that the variability of density with depth should be taken into account when analyzing the wave phenomena in a liquid [5].

The results of the analysis of the first papers, which considered the variability of the density of the liquid in depth, were presented [6], but for some unknown reason, mathematicians and mechanics did not pay attention to these topics. For example, such a fundamental characteristic as the buoyancy frequency calculated in [7], which is the limiting frequency for propagating internal waves in a continuously stratified fluid, escaped from attention. Independently the buoyancy frequency was re-discovered as the natural oscillation frequency of probe balls drifting in a stratified atmosphere [8] and later it was associated with the local extreme in the spectrum of high-frequency pressure oscillations recorded by a microbarograph [9].

Due to the practical importance of the issue, for a long time, researchers have been concentrating their efforts on studying periodic phenomena in a liquid in the context of the force action of waves on obstacles and the wave resistance of bodies moving in a liquid. The scientists assumed the constancy of density and incompressibility of the liquid [10] to be the most important qualities. Nowadays the constant density approximation is still considered high priority in the description of wave processes [11].

In theoretical works much attention is paid to the study of nonlinear properties of waves. The first heuristic discovery—the theory of periodic nonlinear vortex waves [12]—is

**Citation:** Chashechkin, Y.D.; Ochirov, A.A. Periodic Waves and Ligaments on the Surface of a Viscous Exponentially Stratified Fluid in a Uniform Gravity Field. *Axioms* **2022**, *11*, 402. https://doi.org/10.3390/ axioms11080402

Academic Editor: Valery Y. Glizer

Received: 11 July 2022 Accepted: 10 August 2022 Published: 15 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

significant in wave theory development. The approach was extended in [13,14], then analyzed in [15] and it has recently been supplemented with new results [16,17]. Another large cycle of studies of nonlinear waves, initiated by observations of a solitary wave in a shipping channel [18], is successfully proceeding at the present time. The number of various model equations [19,20] for the description of nonlinear wave properties is continuously increasing [21].

Seminal works [22,23] occupied a special place among the first publications. In these studies, the parameters of infinitesimal waves were determined, the waves of finite amplitude were calculated and showed the existence of the wave transfer of matter (Stokes drift) using the methods of the regular perturbations theory. The mass transfer is caused by nonlinear effects, which are distinctly revealed in the deviation of the waveform from the ideal one, i.e., in a nonlinear wave, the crests become sharper, and the troughs become flatter. The limiting angle between the tangents to the right and left sides of the crest is calculated as *α* = 120◦ in [23], and as *α* = 90◦ in [13]. It depends on the nonlinearity parameter ε = ω2*ζ*0/*g* [24] when describing the propagating potential waves of finite amplitude by means of Lambert's complex functions (*ζ*<sup>0</sup> and ω are the amplitude and frequency of the wave, *g* is the gravity acceleration).

Calculations of the viscous attenuation of waves in a deep liquid and a channel of finite depth were carried out by asymptotic methods in a linear and nonlinear formulation [6,25,26]. Such calculations continue to be conducted at the present time using various approximations [27–30].

The "boundary layer" ideas (i.e., a mathematical model of the flow adjacent to the surface of the liquid, in which the influence of viscosity is considerably revealed), which were formulated at the beginning of the 20th century, had a significant impact on the development of the theory of waves in a liquid. The approximation of the boundary layer, accompanied by the reformulation of the defining equations [31], stimulated the development of the search for new analytical methods of finding their solutions [32]. The analysis of the wave equations by the theory of regular perturbations methods showed that in a homogeneous liquid with kinematic viscosity ν, a surface gravitational wave with a frequency <sup>ω</sup> is accompanied by a boundary layer with a specific thickness *<sup>δ</sup>*<sup>ω</sup> <sup>=</sup> <sup>√</sup>2ν/ω. Separate boundary layers are formed on the free surface and on the solid bottom of the channel through which the waves propagate [33].

The boundary layer greatly influences all the parameters of fluid flows, namely pressure distribution, velocity and substance transfer characteristics [34]. Corrections of the theory formulas [33] due to the viscous attenuation of waves were later calculated in the second order of the perturbation theory [35]. The analysis of the dispersion relation for surface gravitational waves and accompanying boundary layers in a viscous liquid, supplemented by calculations of the attenuation coefficient and the scale of spatial attenuation of the wave was carried out in [36].

A more complex model of a "double boundary layer", inside which there is a periodic boundary layer with a length scale *<sup>δ</sup>*<sup>ω</sup> <sup>=</sup> <sup>√</sup>2ν/ω, was proposed to describe standing waves [37].

With the frequency increasing, the surface tension influences more the pattern of waves on the surface of the liquid. The surface tension coefficient *σ* is a parameter that determines the type of dispersion relation for short waves [38] and accompanying fine constituents [39]. The surface tension changes the dependences of the group and phase velocity, the attenuation coefficient and the velocity of matter transfer *U*<sup>ρ</sup> on frequency ω, wave vector **k** or wavelength *λ*. The pattern of waves generated by a short-range localized source in a viscous liquid, taking into account the effects of surface tension, was calculated in [40].

The data of the first systematic experimental studies of surface waves, which were conducted in laboratory basins at the end of the 19th century [41], generally showed satisfactory agreement with the calculations of velocity in the liquid thickness [23]. The visualization of the velocity profile using a "marker" (i.e., a colored wake of a submerging particle) showed the existence of a drift in the near-surface and bottom layers in the direction of wave propagation and slow counterflow in the middle of the liquid layer [42]. Additional cleaning from dust and film of the target liquid surface and consideration of viscous attenuation in the calculations of the near-surface boundary layer significantly reduced the difference between the data of calculations and experiments in a laboratory channel [43]. The observations of rearrangement of the distribution pattern of initially homogeneous suspension in the field of standing gravity waves in a vertically oscillating basin are given in [44].

In the field of propagating gravitational-capillary waves (5 < ω < 50 Hz), the substance of the colored drop is distributed unevenly over the surface of the primary cavity [45], the same as in the case of a drop merging with a target liquid at rest [46]. In the evolutionary process a slowly drifting colored primary contact area, a near-surface jet with a vortex head and a sinking ring vortex remain in the distribution pattern of the substance. There is a pronounced fibrous structure in the distribution of the drop substance in the target fluid [45].

Widespread capillary waves, which are observed in various rivers, seas and oceans, are caused by weak wind gusts in calm [47] and are certainly present on the slopes of large gravitational waves [48] in the stormy open ocean [49] and coastal regions [50]. Capillary waves contribute to the formation of bubbles, drops and foam [51] and participate in the processes of generating sound packets when drops fall [52]. Drop impact sounds form rain noise [53], which is one of the main sources determining the acoustic background of the ocean [54].

The capillary ripples determine the roughness and the real area of the contact surface. The estimation of their influence on the transfer of momentum, energy and matter between the atmosphere and the hydrosphere is of scientific and practical interest. A continuously expanding list of scientific tasks supports the interest in studying the dynamics of capillary waves in a wide range of conditions.

The analysis of the additional dissipation of short gravitational and gravitationalcapillary waves caused by "parasitic" capillary waves on the slopes of longer waves in a laboratory pool are given in [55]. The presence of a surfactant film dampens short capillary waves, which in turn amplifies decimeter waves [56]. The discussion of modern optical methods for studying capillary waves and the data of detailed experiments is carried out in [57]. A review of papers on the nonlinear interaction of coexisting waves of different frequencies within the framework of the theory of "wave turbulence" is presented in [58]. The basic mechanism of energy transfer in weak turbulence theory is validated experimentally in the gravity (four-wave interactions) and capillary (three-wave interactions) regimes. Advanced experiments enable the achievement of full spatiotemporal reconstruction of the wave field in a weakly or strong nonlinear regime, to infer wave statistics as well as waves' nonlinear dispersion relationship, to compare with theory [58]. The study in [59] discusses the influence of attenuation on nonlinear interactions of short waves and surface waves and the properties of "wave turbulence". In all cases, the density of the liquid is assumed to be homogeneous.

However, in real conditions, the influence of differences in atmospheric and hydrosphere temperatures, insolation, radiative heat transfer, variability in the concentration of dissolved and suspended particles and flows ensure the heterogeneity of the density profile throughout the depth of the liquid [60]. The same happens in a thin near-surface layer, where temperature changes rapidly with depth and forms a "cold" or "warm" film [61].

The given study of the properties of surface waves and accompanying fine components in a viscous stably stratified liquid is based on a complete system of fundamental equations, which was firstly collected in [62] and later was considered in [63,64]. In this paper, a reduced system of the fundamental equations in which the density variability is preserved without considering variations in temperature, salinity or pressure that cause it to change is used. This system of equations is analyzed by methods of the singular perturbations

theory [65] taking into account the compatibility condition [66], enabling calculation of both the waves themselves and the family of accompanying fine constituents (i.e., ligaments).

Analysis of the solutions of the fundamental equations system in a continuously stratified fluid has shown that for periodic flows (with a frequency less than the buoyancy frequency), regular solutions describe waves. The wave frequency ω is related to the wave number **k** by an algebraic dispersion relation ω = ω(**k**, **Ak**), which may include the amplitude **A** [66,67].

Singular solutions describe ligaments perturbations—extended in some direction and fine in others—outlining wave beams. The transverse scale of the ligaments *δ*<sup>ν</sup> <sup>ω</sup> is determined by the kinematic viscosity of the liquid ν and the wave frequency ω: *δ*<sup>ν</sup> <sup>ω</sup> <sup>=</sup> <sup>√</sup>ν/<sup>ω</sup> (Stokes scale) or equivalent parameter with the buoyancy frequency *N*: *δ*<sup>ν</sup> *<sup>N</sup>* <sup>=</sup> <sup>√</sup>ν/*N*. The number and positions of the ligaments in space are determined by the geometry of the source, the amplitude and frequency of its oscillations. Calculations of internal wave beams and accompanying ligaments in stratified liquids and gases are given in [68,69].

In this paper, a dispersion equation in a linear approximation for two-dimensional periodic perturbations on the surface of a viscous exponentially stratified liquid is constructed for the first time. The properties of its complete solutions are analyzed. Approximate solutions are obtained by perturbation theory. Regular solutions describe waves in the range from infra-low frequency to gravitational-capillary and capillary waves. The dispersion relations continuously transform into the well-known formulas of the theory of waves in a homogeneous viscous or ideal liquid. Singular solutions characterize ligaments—fine constituents that complement waves. Graphs of the dependence of the wavelengths and ligaments on the frequency, phase and group velocity of the constituents on the wavelength are given. The results can be used to solve problems of generation and propagation of surface waves with physically justified initial and boundary conditions and comparison with a high-resolution experiment.

#### **2. Periodic Flows in a Viscous Exponentially Stratified Fluid**

We consider the propagation of periodic perturbations over the surface of a viscous exponentially stratified fluid in a uniform gravity field. The undisturbed liquid filling the lower half-space is regarded to be incompressible; meanwhile the effects of thermal conductivity and diffusion are neglected.

We perform the study in a Cartesian coordinate system *Oxyz* in which the plane *Oxy* coincides with the equilibrium position of a free surface of the liquid, and the axis *Oz* is directed vertically upwards against the direction of the gravity acceleration **g**. The undisturbed density distribution over depth is exponential ρ0(*z*) = ρ<sup>00</sup> exp(−*z*/Λ). It is characterized by reference density ρ00(*z*0) (i.e., the density value at the equilibrium level *z* = *z*0), as well as the scale Λ = |*d* ln ρ/*dz*| −1 , frequency *N* = ,*g*/Λ and buoyancy period *Tb* = 2*π*/*N*.

Mostly, changes of real fluids density are usually small and produce a small impact on the inertial properties of the flows. Nevertheless, the conservation of terms describing the stratification effects in the governing equations set is important since gravity acceleration is large. In this regard, it is useful to consider three types of medium: stratified fluids when buoyancy scale Λ, frequency *N* and period *Tb* which are included in the list of main parameters; then very weakly stratified fluids, when the scale of buoyancy substantially exceeds the values of other length scales of the problem (so called potentially homogeneous fluid) and actually homogeneous liquid whose density is assumed to be constant in the entire space [70]. Using a weak but variable density helps to save the rank of complete non-linear set and order of the linearized set of governing equations [66,71] and analyzes additional solutions which are lost in approximation of homogeneous fluid.

Depending on the magnitude of the density gradient created by the corresponding temperature or salinity distributions, it is customary to distinguish strongly and weakly stratified fluids, as well as potentially and actually homogeneous fluids. The values of the characteristic physical quantities are given in Table 1 [70].


**Table 1.** Values of characteristic physical quantities.

Strong stratification (*Tb* ∼ 10 s) is created in laboratory installations. Relatively weak stratification (*Tb* ∼ 10 min) is observed in natural conditions. Figure 1 shows the density– depth dependencies for exponential stratification (solid line) and for linear stratification (dotted line) for a strongly stratified fluid *N* = 1 s−1.

**Figure 1.** The density–depth dependencies for exponential stratification (solid line) and for linear stratification (dotted line) for strongly stratified fluid (*N* = 1 s−1, ν = 0.01 St, *σ* = 72 dyn/cm, ρ<sup>00</sup> = 1.0 g/cm3).

To compare the properties of exponentially and linearly stratified media, we analyze the changes in the density gradient with depth. For exponential stratification, the magnitude of the density gradient depends on the vertical coordinate *z*.

$$\frac{d\mathbf{p}}{dz} = -\frac{\mathbf{p}\_{00}}{\Lambda} \mathbf{e}^{-\frac{\mathbf{p}}{\Lambda}}\tag{1}$$

For linear stratification, the density gradient does not depend on the depth:

$$\frac{d\rho}{dz} = -\frac{\rho\_{00}}{\Lambda} \tag{2}$$

For small variations *z* Λ the dependencies of the density gradient are not distinguishable for exponential and linear stratification. The exponent can be represented as a Taylor series expansion over a small parameter *z*/Λ and expression (1) will be written as follows:

$$\frac{d\rho}{dz} = -\frac{\rho\_{00}}{\Lambda} \mathbf{e}^{-\frac{z}{\Lambda}} \simeq -\frac{\rho\_{00}}{\Lambda} \left( 1 - \frac{z}{\Lambda} + \frac{z^2}{2\Lambda^2} \right) \tag{3}$$

The numerical values for the stratification scale for a strongly stratified fluid are <sup>Λ</sup> <sup>=</sup> 9.81 · <sup>10</sup><sup>2</sup> cm, and for a weakly stratified fluid <sup>Λ</sup> <sup>=</sup> 9.81 · 106 cm. At depths of less than 10% of the stratification scale, it can be argued with a high degree of accuracy that the calculations in the model of a linearly stratified fluid and exponentially stratified fluid coincide.

We study below two-dimensional periodic flows of the form *A* = *A*<sup>0</sup> exp *i*(**kx** − ω*t*) with a positive definite frequency ω and a complex wave number **k**, the imaginary part of which characterizes the spatial attenuation of the flow. The disturbances which are homogeneous in the direction of the transverse horizontal coordinate *y* are selected. They include the displacement of the free surface position *ζ*(*x*, *t*). The velocity **u** = *ux***e***<sup>x</sup>* + *uz***e***<sup>z</sup>* with horizontal and vertical velocity components *ux*, *uz* in an incompressible fluid (div**u** = 0) is represented by derivatives of the stream function ψ:

$$
\mu\_x = \partial\_z \psi, \qquad \mu\_z = -\partial\_x \psi \tag{4}
$$

The expression for the density of the liquid ρ = ρ00(*r*0(*z*) + *s*(*x*, *z*, *t*)) replacing the equation of state in [66] includes the initial distribution *r*0(*z*) and the perturbation caused by the examined periodic flow *s*(*x*, *z*, *t*).

Taking into account these assumptions, the system of continuity and Navier–Stokes equations is extremely simplified and takes the traditional form [62,63,69]:

$$z < \mathbb{Q}: \begin{cases} \rho = \rho\_{00}(r(z) + s(x, z, t)) \\ \rho \partial\_t \mathbf{u} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} - \rho \mathbf{v} \Delta \mathbf{u} = \rho \mathbf{g} - \nabla P \\ \partial\_t \rho + \text{div}(\rho \mathbf{u}) = 0 \end{cases} \tag{5}$$

Kinematic and dynamic boundary conditions on a perturbed surface are traditional [62]:

$$z = \mathbb{\tilde{J}} : \begin{cases} \partial\_l(z - \mathbb{\tilde{J}}) + \mathbf{u} \cdot \nabla(z - \mathbb{\tilde{J}}) = 0 \\ \mathbf{\tau} \cdot (\mathbf{n} \cdot \nabla \mathbf{u}) + \mathbf{n} (\mathbf{\tau} \cdot \nabla \mathbf{u}) = 0 \\ P - P\_0 - \sigma \nabla \cdot \mathbf{n} - 2\rho \mathbf{v} \mathbf{n} (\mathbf{n} \cdot \nabla \mathbf{u}) = 0 \end{cases} \tag{6}$$
 
$$\mathbf{n} = \frac{\nabla(z - \mathbb{\tilde{J}})}{|\nabla(z - \mathbb{\tilde{J}})|} = \frac{-\partial\_\tau \mathbb{\tilde{J}} \mathbf{e}\_\mathbf{x} + \mathbf{e}\_\mathbf{z}}{\sqrt{1 + \left(\partial\_\tau \mathbb{J}\right)^2}}, \qquad \qquad \qquad \qquad \qquad \qquad \tau = \frac{\mathbf{e}\_\mathbf{x} + \partial\_\tau \mathbb{\tilde{J}} \mathbf{e}\_\mathbf{z}}{\sqrt{1 + \left(\partial\_\tau \mathbb{J}\right)^2}},$$

where *P* is the hydrodynamic pressure, *σ* = *γuprho*<sup>00</sup> is the total coefficient of the surface tension of the liquid, **n** and τ are the vectors of the external normal and tangent to the free surface of the liquid, respectively.

The studied problem includes the following dimensional parameters: density ρ and its gradient *<sup>d</sup>*<sup>ρ</sup> *dz* , dynamic *μ* and density-normalized kinematic viscosity ν = *μ*/ρ, gravity acceleration *g*, surface tension coefficient *σ* and density-normalized surface tension coefficient *γ* = *σ*/ρ, amplitude *A*, frequency ω period *T*<sup>ω</sup> = 2*π*/ω and wavelength *λ*. The wave is also characterized by a wave vector **k**, the modulus of which is related to the wavelength *λ* = 2*π*/|**k**|.

In the transition to a mathematical description, physical quantities serve as the basis for the introduction of characteristic scales of length, time and speed. These scales are used in the physical interpretation of the results obtained, assessing the degree of influence of various physical factors.

The degree of relative influence of viscosity and gravity on fluid flows is characterized by the scale *δ*<sup>ν</sup> *<sup>g</sup>* <sup>=</sup> ,<sup>3</sup> <sup>ν</sup>2/*<sup>g</sup>* that appears in [71].

Traditionally, stratification is characterized by its own length scale which takes into account or neglect the effects of compressibility, the frequency *<sup>N</sup>* = ,*g*/Λ, period *Tb* = <sup>2</sup>*π*/*<sup>N</sup>* and the inverse value of the buoyancy frequency *TN* = *Tb*/2*π*. The stratification itself is characterized by a combinational viscous wave scale *Lg*<sup>ν</sup> *<sup>N</sup>* = (*g*ν) 1/3/*N* [72,73], as well as a length scale *δ*<sup>ν</sup> *<sup>N</sup>* <sup>=</sup> <sup>√</sup>ν/*N*– a functional analogue of the dissipative Stokes microscale *δ*ν <sup>ω</sup> <sup>=</sup> <sup>√</sup>ν/ω.

The dualism of the nature of surface waves, whose properties depend on the gravity acceleration *g*, normalized coefficients of surface tension *γ* and kinematic viscosity ν results in the introduction of several scales, including the capillary length *δ γ <sup>g</sup>* = ,*γ*/*g*. The inequality *λ* < *δ γ <sup>g</sup>* indicates the severity of surface force effects, and *λ* >> *δ γ <sup>g</sup>*—indicates the severity of gravitational effects.

The structure of near-surface flows is also characterized by a known proper velocity scale v*<sup>c</sup>* = *<sup>γ</sup>*/<sup>ν</sup> and an additional time scale of kinematic nature *<sup>T</sup><sup>γ</sup>* <sup>ν</sup>*<sup>g</sup>* = v*c*/*g* = *γ*/ν*g*.

Wave propagation is characterized by group **c***<sup>g</sup>* = *<sup>∂</sup>*<sup>ω</sup> *<sup>∂</sup>***<sup>k</sup>** and phase velocity **<sup>c</sup>***ph* <sup>=</sup> <sup>ω</sup>**<sup>k</sup> k**2 . Ratios of uniform scales form a set of dimensionless parameters, parts of which are

used in further calculations.

#### *2.1. Equations of Periodic Flows and Dispersion Relations for Plane Infinitesimal Waves*

To simplify expressions, we accept the initial density distribution to be exponential, as it is customary in most models:

$$
\rho\_0(z) = \rho\_{00} r(z) = \rho\_{00} \exp(-z/\Lambda) \tag{7}
$$

The change of variables proposed in [74] allows the transformation of the equations of motion for an arbitrary smooth density profile into equations with constant coefficients that determine the preferential choice of the exponential density profile.

Further calculations are carried out in the Boussinesq approximation, considering the smallness of absolute density variations in stratified media (in particular, a small value of the wavelength and buoyancy scale ratio C = *λ*/Λ). The density variations are neglected everywhere, except for the term with a large coefficient (i.e., the gravity acceleration **g**). In this approximation, the system (5) takes the form:

$$\begin{cases} \partial\_t \boldsymbol{\Psi} + \partial\_x \boldsymbol{\Psi} \partial\_x \boldsymbol{\Psi} - \partial\_x \boldsymbol{\Psi} \partial\_{zz} \boldsymbol{\Psi} - \nu \partial\_x \boldsymbol{\Psi} \boldsymbol{\Psi} + \partial\_x P = 0 \\\ \boldsymbol{\mathcal{g}}(\boldsymbol{r} + \mathbf{s}) - \boldsymbol{\mathcal{g}} - \partial\_{tx} \boldsymbol{\Psi} + \partial\_x \boldsymbol{\Psi} \partial\_{xz} \boldsymbol{\Psi} - \partial\_z \boldsymbol{\Psi} \partial\_{xx} \boldsymbol{\Psi} - \nu \partial\_x \boldsymbol{\Lambda} \boldsymbol{\Psi} + \partial\_z P = 0 \\\ \partial\_t \boldsymbol{s} + \partial\_z \boldsymbol{\Psi} \partial\_x \boldsymbol{s} - \partial\_x \boldsymbol{\Psi} \partial\_z (\boldsymbol{r} + \mathbf{s}) = 0 \end{cases} \tag{8}$$

In a linear approximation the system (8) reduces to:

$$\begin{cases} \partial\_{tz}\psi - \mathbf{v}\partial\_z\Delta\psi + \partial\_x P = 0\\ \mathbf{g}(r+s) - \mathbf{g} - \partial\_{tx}\psi - \mathbf{v}\partial\_x\Delta\psi + \partial\_z P = 0\\ \partial\_t s - \partial\_x \psi \partial\_z r = 0 \end{cases} \tag{9}$$

Cross-differentiation of spatial coordinates of the upper and middle equations of the system (9) allows getting rid of pressure:

$$\begin{cases} \left. \frac{\partial\_t \Delta \psi}{\partial t} - \mathcal{g} \partial\_x \mathbf{s} - \mathbf{v} \Delta \Delta \psi = 0 \\\ \partial\_t \mathbf{s} - \partial\_x \psi \partial\_z r = 0 \end{cases} \tag{10}$$

Subtraction of the second equation, differentiated by coordinate and multiplied by a coefficient *g* from the first equation of the system (10), differentiated by time, allows getting rid of the function *s*(*x*, *z*, *t*):

$$
\partial\_{tt} \Delta \psi - \mathbf{v} \partial\_t \Delta \Delta \psi + \partial\_{xx} \psi \partial\_z r = 0 \tag{11}
$$

Equation (11) takes the next form for an exponentially stratified fluid (7):

$$
\partial\_{tt} \Delta \psi - \mathbf{v} \partial\_t \Delta \Lambda \psi + N^2 \exp(-z/\Lambda) \partial\_{xx} \psi = 0, \quad N = \sqrt{\mathbf{g}/\Lambda} \tag{12}
$$

The boundary conditions (6) for all variables of the infinitesimal waves *A* = *A*<sup>0</sup> exp(*i***kx** − *i*ω*t*) are traditionally carried away from the wave surface *z* = *ζ* to the unperturbed level *z* = 0 and take the form

$$z = 0 \\ \begin{cases} \partial\_t \zeta + \partial\_x \psi = 0 \\ -\rho g \zeta + \rho\_{00} P + 2\rho\_{00} \mathbf{v} \partial\_{xx} \psi + \sigma \partial\_{xx} \zeta = 0 \\ \partial\_{zz} \psi - \partial\_{xx} \psi = 0 \end{cases} \tag{13}$$

The small coefficient (kinematic viscosity ν), which ensures slow attenuation of the wave as it propagates in the liquid, allows the classification of the system (12) as a system of singularly perturbed equations and the application of the theory of singular perturbations [65] for its analysis. The compatibility condition, which requires the analysis of all the roots of the dispersion equation and the system (12) solutions, should be taken into account [66].

For surface waves with ω > *N*, the solution of Equation (12) is sought in the form of plane waves:

$$
\Psi = \left( A\_{+} e^{ik\_{x}x - i\omega t} + A\_{-} e^{-ik\_{x}x - i\omega t} \right) \left( e^{k\_{z}z} + \beta e^{k\_{l}z} \right) \tag{14}
$$

Here *kz* corresponds to the regular solution of the dispersion equation, and *kl* corresponds to the singular solution of the dispersion equation.

We obtain the shape of the deviation of the free surface of the liquid, the condition of the relationship between the amplitudes of the regular and singular component by substituting (14) into Equation (12) and boundary conditions (13):

$$\begin{cases} \mathcal{J} = \frac{k\_x}{\omega \nu} \left( A\_+ e^{i\mathbf{k}\_x \mathbf{x} - i\omega \nu t} - A\_- e^{-i\mathbf{k}\_x \mathbf{x} - i\omega \nu t} \right) (1 + \beta) \\ \mathcal{J} = -\frac{k\_x^2 + k\_z^2}{k\_x^2 + k\_l^2} \end{cases} \tag{15}$$

as well as the dispersion relations:

$$\begin{cases} -N^2 k\_x^2 + e^{z/\Lambda} \left( k\_x^2 - k\_{z,l}^2 \right) \omega \left( ik\_x^2 \mathbf{v} - ik\_{z,l}^2 \mathbf{v} + \omega \right) = 0\\\ g k\_x^2 - 3ik\_x^2 k\_{z,l} \mathbf{v} \omega + k\_{z,l} \left( ik\_{z,l}^2 \mathbf{v} - \omega \right) \omega + k\_x^4 \gamma = 0, \qquad \gamma = \sigma/\rho\_{00} \end{cases} \tag{16}$$

Calculated for the first time the expressions (16) are transformed into the relations for a viscous homogeneous fluid and for an ideal exponentially stratified fluid in limiting transitions *N* → 0 and ν → 0 respectively. Similar dispersion relations for internal waves and ligaments in the thickness of a stratified fluid were presented in [75].

#### *2.2. Solution of the Dispersion Equation*

It is convenient to analyze the dispersion relations (13) in a dimensionless form. The time scale will be the parameter, which is inverse to the buoyancy frequency *N*, and the viscous wave scale *L*<sup>ν</sup> = (ν*g*) 1/3/*N* is chosen as the length scale. The ratios of the eigenscales of the problem *δ* = *<sup>δ</sup> γ g δ*ν *N* <sup>2</sup> = *<sup>γ</sup> <sup>g</sup>* · *<sup>N</sup> <sup>v</sup>* and <sup>ε</sup> <sup>=</sup> *Lv* <sup>Λ</sup> <sup>=</sup> *<sup>N</sup>*ν1/3 *<sup>g</sup>*2/3 are used to construct new dimensionless parameters *δ*, ε involved in further calculations. Then, expressions (16) could be rewritten in a dimensionless form:

$$\begin{cases} -k\_{\ast x}^2 + i\varepsilon^{z/\Lambda} \left(k\_{\ast x}^2 - k\_{\ast z,l}^2\right)^2 \varepsilon \boldsymbol{\omega}\_{\ast} + \varepsilon^{z/\Lambda} \left(k\_{\ast x}^2 - k\_{\ast z,l}^2\right) \boldsymbol{\omega}\_{\ast}^2 = 0\\\ k\_{\ast x}^2 + k\_{\ast x}^4 \delta \varepsilon - 3ik\_{\ast x}^2 k\_{\ast z,l} \varepsilon^2 \boldsymbol{\omega}\_{\ast} + ik\_{\ast z,l}^3 \varepsilon^2 \boldsymbol{\omega}\_{\ast} - k\_{\ast z,l} \varepsilon \boldsymbol{\omega}\_{\ast}^2 = 0\\\ \delta = N\gamma/\gamma \text{g}\_{\prime} & \varepsilon = N\mathbf{v}^{1/3}/\text{g}^{2/3} \end{cases} \tag{17}$$

The upper dispersion relation in (17) has regular solutions, which we have denoted *k*∗*z*, and singular solutions, which are denoted *k*∗*l*. For a large number of real liquids, the parameter ε 1 and regular solutions are found by direct decomposition over a small parameter ε:

$$k\_{\ast z} = k\_{\ast 0z} + \varepsilon k\_{\ast 1z} + \varepsilon^2 k\_{\ast 2z} + \dots \tag{18}$$

Substituting (18) into the upper expression in (17) we obtain up to the terms of the order *O* # ε2 \$ :

$$k\_{\ast z} = \pm \frac{\sqrt{\omega\_{\ast}^2 - e^{-z/\Lambda}}}{\omega\_{\ast}} k\_{\ast x} \pm \varepsilon \frac{i e^{-2z/\Lambda} k\_{\ast x}^3}{2 \omega\_{\ast}^4 \sqrt{\omega\_{\ast}^2 - e^{-z/\Lambda}}} \tag{19}$$

Since the equations in the system (17) have the fourth degree, it is necessary to find two more solutions, which are singular perturbed ones in the form of decomposition [65]:

$$k\_{\ast I} = \varepsilon^{-\eta} \left( k\_{\ast 0I} + \varepsilon k\_{\ast 1I} + \varepsilon^2 k\_{\ast 2I} + \dots \right) \tag{20}$$

The parameter η in (20) is chosen in such a way that at the highest degree the main term of the decomposition persists. Substituting instead of *k*∗*<sup>l</sup>* in the upper expression in (17) *<sup>k</sup>*∗*<sup>l</sup>* <sup>=</sup> <sup>ε</sup>−η*k*∗0*<sup>l</sup>* and equating the exponents <sup>ε</sup> in different terms we get:

$$\begin{array}{l} 1 - 4\eta = -2\eta \\ 1 - 4\eta = 1 - 2\eta \\ 1 - 4\eta = 0 \\ 1 - 2\eta = -2\eta \\ 1 - 2\eta = 0 \\ -2\eta = 0 \end{array} \tag{21}$$

At the highest degree, the main term of the decomposition remains only at η = 1/2. Substituting the value η = 1/2 in (20) and then in (17) up to the terms of the order *O* ε3/2 we get:

$$k\_{\ast l} = \pm \frac{(1 - i)\sqrt{\omega\_{\ast}}}{\sqrt{2\varepsilon}} \pm \frac{(1 + i)\left(\varepsilon^{-z/\Lambda} + \omega\_{\ast}^{2}\right)k\_{\ast \ast}^{2}}{2\sqrt{2}\omega\_{\ast}^{5/2}}\sqrt{\varepsilon} \tag{22}$$

Leaving only the main terms of ε in the lower dispersion relation (17), we obtain the dispersion equation for the regular part of the solution:

$$k\_{\ast x} + \left(k\_{\ast x}^3 \delta - \omega\_{\ast}\sqrt{\omega\_{\ast}^2 - 1}\right)\varepsilon = 0\tag{23}$$

For the singular part of the solution, leaving only the main terms of the expansion of ε, we obtain:

$$1 + k\_{\ast x}^2 \delta \varepsilon + \frac{(1+i)\left(1 - 2\omega\_{\ast}^2\right)}{\sqrt{2\omega\_{\ast}}} \varepsilon^{3/2} = 0\tag{24}$$

The solution of the dispersion Equation (23) for the regular wave part is:

$$\begin{aligned} k\_{\ast x} &= -\frac{2^{1/3}}{\mathfrak{a}} + \frac{\mathfrak{a}}{3 \cdot 2^{1/3} \delta \varepsilon} \\ k\_{\ast x} &= \frac{1 \pm i\sqrt{3}}{2^{2/3} \mathfrak{a}} - \frac{\left(\frac{1 \mp i\sqrt{3}}{6 \cdot 2^{1/3} \delta \varepsilon}\right) \mathfrak{a}}{6 \cdot 2^{1/3} \delta \varepsilon} \\ \mathfrak{a} &= \left(27 \delta^2 \varepsilon^3 \mathfrak{w}\_\* \sqrt{\mathfrak{w}\_\*^2 - 1} + \sqrt{108 \delta^3 \varepsilon^3 + 729 \delta^4 \varepsilon^6 \mathfrak{w}\_\*^2 (\mathfrak{w}\_\*^2 - 1)}\right)^{1/3} \end{aligned} \tag{25}$$

From the condition of the physical realization of the solution (i.e., the damping of the flow with depth) it follows that only roots with Re(*k*∗*z*) > 0 possess a physical meaning. Decomposing the solution (25) for the regular wave part into a Taylor series, we obtain:

$$\begin{aligned} k\_{\ast x} &= \omega\_{\ast} \varepsilon \sqrt{\omega\_{\ast}^{2} - 1} + O(\varepsilon^{4}) \\ k\_{\ast x} &= \pm \frac{i}{\sqrt{\delta \varepsilon}} - \omega\_{\ast} \varepsilon \frac{\sqrt{\omega\_{\ast}^{2} - 1}}{2} + O\left(\varepsilon^{5/2}\right) \end{aligned} \tag{26}$$

With consideration of (19) and the condition Re(*k*∗*z*) > 0, it can be seen that we implement only one root of (26) and finally for the regular wave solution we get:

$$k\_{\ast X} = -\frac{2^{1/3}}{\mathfrak{a}} + \frac{\mathfrak{g}}{3 \cdot 2^{1/3} \delta \varepsilon}$$

$$k\_{\ast z} = \frac{\sqrt{\omega\_\*^2 - \varepsilon^{-z/\Lambda}}}{\omega\_\*} \left( -\frac{2^{1/3}}{\mathfrak{a}} + \frac{\mathfrak{g}}{3 \cdot 2^{1/3} \delta \varepsilon} \right) + \varepsilon \frac{\mathfrak{a}^{-2z/\Lambda} \left( -\frac{2^{1/3}}{\mathfrak{a}} + \frac{\mathfrak{a}}{3 \cdot 2^{1/3} \delta \varepsilon} \right)^3}{2 \omega\_\*^4 \sqrt{\omega\_\*^2 - \varepsilon^{-z/\Lambda}}} \tag{27}$$

We perform similar transformations for singular ligament-solution. The solutions (24), which are found taking into account (22) and the conditions of the physical implementation of the roots Re(*k*∗*l*) > 0 of the singular ligament, take the form:

$$\begin{array}{l}k\_{\ast x} = \pm \sqrt{\left(\delta \varepsilon\right)^{-1} \left(\left((1+i)(2\omega\_{\ast}^{2}-1)/\sqrt{2\omega\_{\ast}}\right) \varepsilon^{3/2} - 1\right)}\\k\_{\ast I} = \frac{(1-i)2\delta\omega\_{\ast}^{3} + (1+i)\left(\varepsilon^{-z/\Lambda} + \omega\_{\ast}^{2}\right)\left(\left((1+i)\left(2\omega\_{\ast}^{2}-1\right)/\sqrt{2\omega\_{\ast}}\right) \varepsilon^{3/2} - 1\right)}{2\sqrt{2}\omega\_{\ast}^{5/2}\delta\sqrt{\varepsilon}}\end{array} \tag{28}$$

Expressions (27), (28) describe new dispersion relations for surface waves (27) and associated ligaments (28) in a viscous exponentially stratified fluid.

#### *2.3. Low Frequency Waves*

The solution for low frequency waves (ω < *N*) is sought as the sum of the propagation of gravitational waves (waves with wave number components *kz*) and fine perturbation (with wave number component *ks*) as well:

$$\Psi = \left( A\_{+} e^{i\mathbf{k}\_{x}\mathbf{x} - i\omega\mathbf{v}t} + A\_{-} e^{-i\mathbf{k}\_{x}\mathbf{x} - i\omega\mathbf{v}t} \right) \left( a e^{i\mathbf{k}\_{z}z} + \beta e^{-i\mathbf{k}\_{z}z} + \chi e^{\mathbf{k}\_{z}z} \right) \tag{29}$$

Substituting (26) into (9) leads to dispersion relations:

$$\begin{cases} -N^2 k\_x^2 + \varepsilon^{z/\Lambda} \left(k\_x^2 + k\_z^2\right) \omega \left(ik\_x^2 \mathbf{v} + ik\_z^2 \mathbf{v} + \omega\right) = 0\\ -N^2 k\_x^2 + i\varepsilon^{z/\Lambda} \left(k\_s^2 - k\_x^2\right) \omega \left(k\_s^2 \mathbf{v} - k\_x^2 \mathbf{v} + i\omega\right) = 0 \end{cases} \tag{30}$$

Substituting solution (29) into the boundary conditions (13), we obtain an expression for the shape of the free surface and additional relations. These additional relations determine the relationship between the amplitudes of the various components of periodic motion:

$$\mathcal{J} = \left( A\_{+} e^{ik\_{x}x - i\omega\nu t} - A\_{-} e^{-ik\_{x}x - i\omega\nu t} \right) (a + \beta + \chi) k\_{x}/\omega \tag{31}$$

$$(\alpha + \beta) \left(k\_x^2 - k\_z^2\right) + \chi \left(k\_x^2 + k\_s^2\right) = 0\tag{32}$$

$$\lambda k\_x^2 \left(\xi + \gamma k\_x^2\right) \left(\mathfrak{a} + \mathfrak{B} + \chi\right) + \text{var}\left(k\_z \mathfrak{a} - \mathfrak{B}\right) \left(3k\_x^2 + k\_z^2\right) + i\chi k\_s \left(k\_s^2 - 3k\_x^2\right)\right) - \mathfrak{w}^2 \left(i k\_z \mathfrak{a} - \mathfrak{B}\right) - \chi k\_s \mathfrak{b} = 0 \tag{33}$$

In a dimensionless form, the dispersion relations (30) take the form:

$$\begin{cases} \begin{aligned} -k\_{\ast x}^{2} + i\varepsilon^{z/\Lambda} \left(k\_{\ast x}^{2} - k\_{\ast s}^{2}\right)^{2} \varepsilon \boldsymbol{\omega}\_{\ast} + \varepsilon^{z/\Lambda} \left(k\_{\ast x}^{2} - k\_{\ast s}^{2}\right) \boldsymbol{\omega}\_{\ast}^{2} = 0 \\\ -k\_{\ast x}^{2} + i\varepsilon^{z/\Lambda} \left(k\_{\ast x}^{2} + k\_{\ast z,l}^{2}\right)^{2} \varepsilon \boldsymbol{\omega}\_{\ast} + \varepsilon^{z/\Lambda} \left(k\_{\ast x}^{2} + k\_{\ast z,l}^{2}\right) \boldsymbol{\omega}\_{\ast}^{2} = 0 \end{aligned} \tag{34}$$

The expressions for the surface component of periodic motion, up to the notation, coincide with the dispersion equations of surface waves discussed above. For low frequency waves, we will look for regular wave solutions in the form of expansion (18) and up to the terms of the order *O* # ε2 \$ :

$$k\_{\ast z} = \pm \frac{\sqrt{\varepsilon^{-z/\Lambda} - \omega\_{\ast}^2}}{\omega\_{\ast}} k\_{\ast x} \mp \varepsilon \frac{i e^{-2z/\Lambda} k\_{\ast x}^3}{2 \omega\_{\ast}^4 \sqrt{\omega\_{\ast}^2 - \varepsilon^{-z/\Lambda}}} \tag{35}$$

Similarly to surface waves case, we find singular solutions of the dispersion equation for low frequency waves using decomposition (20). Just as in the case of high frequency surface waves ω > *N*, the exponent of the degree η = 1/2 is the only one that satisfies the condition of the main term of the decomposition at the highest degree *kl*. Making a substitution up to the terms of the order *O* ε3/2 we get:

$$k\_{\ast l} = \pm \frac{(1+i)\sqrt{\omega\_{\ast}}}{\sqrt{2\varepsilon}} \mp \frac{(1-i)\left(\varepsilon^{-z/\Lambda} + \omega\_{\ast}^{2}\right)k\_{\ast x}^{2}}{2\sqrt{2}\omega\_{\ast}^{5/2}}\sqrt{\varepsilon} \tag{36}$$

From the relation (32) we obtain that

$$\chi = (\mathfrak{a} + \mathfrak{z}) \frac{k\_z^2 - k\_x^2}{k\_x^2 + k\_s^2} \tag{37}$$

$$
\alpha + \beta + \chi = (\alpha + \beta) \frac{k\_s^2 + k\_z^2}{k\_x^2 + k\_s^2} \tag{38}
$$

From the ratio (33) follows:

$$\mathbf{a} + \boldsymbol{\mathfrak{B}} = -\frac{2i(k\_s^2 + k\_x^2)kza\left(3k\_x^2\mathbf{v} + k\_z^2\mathbf{v} - i\boldsymbol{w}\right)\boldsymbol{w}}{ik\_x^2(k\_s^2 + k\_z^2)(\boldsymbol{g} + k\_x^2\boldsymbol{\gamma}) - (k\_s + ik\_z)(-k\_s^2k\_x^2 + 3k\_x^4 + 4ik\_sk\_x^2k\_z + (k\_s^2 + k\_x^2)k\_z^2)\boldsymbol{\gamma}\boldsymbol{w} - (k\_s - ik\_z)(ik\_x^2 + k\_zk\_z)\boldsymbol{w}^2} \tag{39}$$

In the low viscosity approximation, we obtain that

$$
\kappa + \beta \simeq 0, \quad \chi \simeq -(\kappa + \beta) = 0 \tag{40}
$$

The relations (40) correspond to the situation when all the energy is concentrated in gravitational waves.

#### *2.4. Periodic Flows on the Surface of a Viscous Exponentially Stratified Liquid*

First of all, let us consider the dependence of the wavelength *λ* on the frequency of the wave motion ω. We define the wavelength as follows:

$$
\lambda = 2\pi / \sqrt{\mathrm{Re}(k\_x^2) + \mathrm{Im}(k\_z^2)}\tag{41}
$$

This method of determination of the wavelength is due to the fact that the imaginary part of the wave number *kx* component and the real part of the wave number *kz* component are responsible for the spatial attenuation of motion and are not impactful to wave propagation. Substituting the dispersion relations (27) in (32), we can construct the desired dependencies. Figure 2 shows the dependence of the wavelength on the frequency of wave motion *λ*(ω).

The velocity of movement of the phase front of structures (wave and ligament) and the rate of energy transfer are of particular interest. The phase front moves with the phase velocity of the wave (and its analogue for the singular solution):

$$\mathbf{c}\_{plt} = \omega \mathbf{k} / |k|^2 \tag{42}$$

The energy is transferred with a group velocity *cgr*, which is defined as follows:

$$\mathbf{c}\_{\mathcal{S}^r} = \left\{ \left( \frac{\partial k\_x}{\partial \omega} \right)^{-1}, \left( \frac{\partial k\_{z,l}}{\partial \omega} \right)^{-1} \right\} \tag{43}$$

Figure 3 shows the dependencies for the phase (dashed line) and group (solid line) velocities.

**Figure 2.** The dependence of the wavelength on the frequency of wave motion in a viscous exponentially stratified fluid for a periodic solution. The curves indicated by the letter (*W*) are constructed for a liquid with ν = 0.01 St, *σ* = 72 dyn/cm, ρ<sup>00</sup> = 1 g/cm3, and by the letter (*Gl*)—for a liquid with glycerin parameters (ν = 11.746 St, ρ<sup>00</sup> = 1.26 g/cm3, *σ* = 64.7 dyn/cm). The numbers indicate a different degree of stratification. Index (1) corresponds to a weak pycnocline *N* = 0.001 s−1, index (2) to a weakly stratified fluid *N* = 0.01 s−1, and index (3) to a strongly stratified fluid *N* = 1 s−<sup>1</sup> for water and glycerin.

Figures 2 and 3 show that viscosity has a noticeable effect on capillary waves with a wavelength *λ* < *δ γ <sup>g</sup>* , and the stratification influences the waves with frequencies close to the buoyancy frequency *N*. With the advent of stratification, a forbidden part of the spectrum appears for surface waves (with frequencies lower than the buoyancy frequency). With the tending to the stratification frequency, the group velocity goes to zero, and the phase velocity tends to infinity. A similar pattern is observed when electromagnetic waves propagate in waveguides. The size of the waveguide sets a certain critical size, beyond which the electromagnetic wave does not propagate in the waveguide channel. In stratified fluids, this mechanism is not due to external geometric constraints, but is determined by the characteristics of the medium (stratification).

**Figure 3.** *Cont*.

**Figure 3.** Dependences of the phase (dashed lines) and group (solid lines) velocities on the frequency of wave motion (**a**) and on the wavelength (**b**) in a viscous exponentially stratified fluid for a periodic solution. The curves indicated by the letter (*W*) are constructed for a liquid with water parameters, and by the letter (*Gl*)—for a liquid with glycerin parameters. The numbers indicate a different degree of stratification. Indexes (1) and (6) correspond to a weak pycnocline *N* = 0.001 s−1, indexes (2) and (5) correspond to a weakly stratified liquid *N* = 0.01 s−<sup>1</sup> and indexes (3) and (4) correspond to a strongly stratified liquid *N* = 1 s−1.

#### **3. Reduction to Approximation of Actually Homogeneous Fluid**

The fine constituents of periodic flows—ligaments—are a consequence of the dissipative properties of the medium, which exist not only in stratified, but also in homogeneous liquids. They are described by singular solutions of dispersion equations, the appearance of which can be observed experimentally in the structure of flows in an inhomogeneous medium. If the effects associated with buoyancy are neglected, then the density of the liquid can be considered as actually homogeneous:

$$
\mathfrak{p} = \mathfrak{p}\_{00} \equiv \text{const} \tag{44}
$$

and the basic equations of motion are simplified. For a viscous homogeneous liquid in a 2D formulation in a linear approximation equations for the stream function are shortened to:

$$
\partial\_t \Delta \psi - \mathbf{v} \Delta \Delta \psi = 0 \tag{45}
$$

Equation (45) can also be obtained by performing a limiting transition *N* → 0 in (12). The boundary conditions that are removed to the equilibrium surface in the case of a homogeneous liquid are as follows:

$$z = 0 \left\{ \begin{array}{l} \partial\_t \zeta + \partial\_x \psi = 0 \\ -g\zeta + P + 2\nu \partial\_{xx} \psi + \gamma \partial\_{xx} \zeta = 0 \\ \partial\_{zz} \psi - \partial\_{xx} \psi = 0 \end{array} \right. \tag{46}$$

Similarly to the basic equations of motion, the boundary conditions (46) can also be obtained from the boundary conditions (13) using a limit transition *N* → 0. Substitution of the solution in the form of a propagating wave

$$\psi = A\_{+} \exp(k\_{z}z + ik\_{x}\mathbf{x} - i\omega t) + A\_{-} \exp(k\_{z}z - ik\_{x}\mathbf{x} - i\omega t) \tag{47}$$

in (45) leads to dispersion relations between the components of the wave number:

$$
\left(k\_x^2 - k\_z^2\right)\left(\mathbf{v}\left(k\_x^2 - k\_z^2\right) - i\omega\right) = 0\tag{48}
$$

The relation (48) naturally decomposes into two solutions:

$$\begin{cases} k\_z^2 = k\_x^2\\ k\_l^2 = k\_x^2 + \frac{i\omega}{\gamma} \end{cases} \tag{49}$$

Here, the reassignment of one solution *kz* into *kl* is introduced. The component *kz* corresponds to the wave solution, and the component *kl* corresponds to the ligament solution. Taking into account (49), the solution of the problem in the form of a propagating wave is transformed into:

$$\psi = \left( A\_{+} \exp(i k\_{x} \mathbf{x} - i \omega t) + A\_{-} \exp(-i k\_{x} \mathbf{x} - i \omega t) \right) \left( \exp(k\_{z} z) + B \exp(k\_{l} z) \right) \tag{50}$$

Substituting (50) into the kinematic boundary condition, we obtain the relationship between the amplitudes of the velocity field and the deviation of the free surface from the equilibrium value:

$$\mathcal{J} = \frac{k\_{\text{x}}}{\omega} (1 + B)(A\_{\text{-}} \exp(-ik\_{\text{x}}\mathbf{x} - i\omega\mathbf{t}) - A\_{\text{+}} \exp(ik\_{\text{x}}\mathbf{x} - i\omega\mathbf{t})) \tag{51}$$

and substitution of (50) into the dynamic condition for tangential tensions leads to the expression for the amplitude of the ligament component:

$$B = -\frac{\left(k\_x^2 + k\_z^2\right)}{\left(k\_x^2 + k\_l^2\right)} = -2k\_x^2 \left(2k\_x^2 + \frac{i\omega}{\mathbf{v}}\right)^{-1} \tag{52}$$

Getting rid of the pressure in the dynamic boundary condition and taking normal components, we rewrite it as:

$$
\Delta\bar{\boldsymbol{\sigma}}\_{tz}\Delta\boldsymbol{\Psi} - \bar{\boldsymbol{\sigma}}\_{ttz}\boldsymbol{\Psi} + \bar{\boldsymbol{\varrho}}\partial\_{xx}\boldsymbol{\Psi} + 2\boldsymbol{\nu}\partial\_{txxx}\boldsymbol{\Psi} - \boldsymbol{\gamma}\bar{\partial}\_{xxxx}\boldsymbol{\Psi} = \boldsymbol{0} \tag{53}
$$

Substituting the solution (50) into the boundary conditions (53) taking into account (49) leads to the dispersion relation for the wave constituent:

$$
\gamma k\_z^3 - 2i\mathbf{v}\omega k\_z^2 + gk\_z - \omega^2 = 0\tag{54}
$$

and ligament constituent:

$$\gamma \left( k\_l^2 - \frac{i\omega}{\mathbf{v}} \right)^2 - 2i\mathbf{v}\omega \left( k\_l^2 - \frac{i\omega}{\mathbf{v}} \right) k\_l + \mathbf{g} \left( k\_l^2 - \frac{i\omega}{\mathbf{v}} \right) + \omega k\_l \left( i\mathbf{v} k\_l^2 - \omega \right) = 0 \tag{55}$$

It can be noticed that Equations (54) and (55) are also obtained from (16) with a limit transition *N* → 0. Further we will transfer the dispersion relations (54) and (55) to dimensionless forms in the same way as it has been done in the previous paragraph. We will choose our own viscous scale *δ*<sup>ν</sup> *<sup>g</sup>* <sup>=</sup> ,<sup>3</sup> <sup>ν</sup>2/*<sup>g</sup>* as the length scale, and the ratio as the time scale will be *T<sup>γ</sup>* <sup>ν</sup>*<sup>g</sup>* = *<sup>γ</sup>*/ν*<sup>g</sup>* and the decomposition parameter <sup>ε</sup> = <sup>3</sup> . *δ*ν *g* 6 / # *δ γ g* \$6 = ,<sup>3</sup> (ν4*g*)/*γ*3.

In a dimensionless form, the conventional dispersion Equation [15] for the wave component of the solution will be rewritten taking into account (54):

$$k\_{\ast x}^3 - 2i\varepsilon^2 \omega\_\ast k\_{\ast x}^2 + \varepsilon k\_{\ast x} - \varepsilon^3 \omega\_\ast^2 = 0\tag{56}$$

The solution of Equation (56) is obtained by the standard method and takes the form:

$$k\_x^\* = \frac{2}{3} \text{i}\varepsilon^2 \omega^\* - \frac{\varepsilon \left(3 + 4\varepsilon^3 \omega^{\*2}\right)}{3\beta} + \frac{1}{3}\beta\tag{57}$$

β = 

$$k\_x = \frac{2}{3} \text{i}\varepsilon^2 \text{w}^\* \pm \frac{1}{6} \text{i}\left(\pm \text{i} + \sqrt{3}\right) \mathfrak{H} + \frac{\left(1 \pm i\sqrt{3}\right) \varepsilon\left(3 + 4\varepsilon^3 \text{w}^{\*2}\right)}{6 \cdot 2^{2/3} \mathfrak{H}}$$

$$-8 \text{i}\varepsilon^6 \text{w}^{\*3} + \frac{9}{2} \varepsilon^3 \text{w}^\*(-2\text{i} + 3\text{w}^\*) + \frac{3}{2} \sqrt{3} \sqrt{\varepsilon^3(4 - 32 \text{i} \text{}^6 \text{w}^{\*5} + \varepsilon^3 \text{w}^{\*2}(4 + 9 \text{w}^\*(-4\text{i} + 3 \text{w}^\*))}\right)^{1/3} \tag{58}$$

Decomposition (57) into a Taylor series by a small parameter ε at least up to the summands *O* ε9/2 gives the expression:

$$k\_{\ast z} = \varepsilon^2 \mathfrak{w}\_{\ast}^2 + O\left(\varepsilon^{9/2}\right) \tag{59}$$

Decomposition of the roots (58) into a Taylor series by a small parameter ε at least up to the summands *O* ε9/2 gives the expression:

$$k\_{\ast z} = \pm \frac{1}{8} \mathrm{i} \sqrt{\varepsilon} \left( 8 \pm 4 \varepsilon^{3/2} \omega\_{\ast} (2 + i \omega\_{\ast}) + \varepsilon^3 \omega\_{\ast}^2 \left( 4 - 4i \omega\_{\ast} + 3 \omega\_{\ast}^2 \right) \right) + O \left( \varepsilon^{9/2} \right) \tag{60}$$

As it follows from the condition of the physical realization of the roots Re(*kz*) > 0, considering the positive definite frequency of the wave motion, only one root, which is given by the expression (57), can be physically realized This solution describes the wave part of a periodic motion in a liquid.

To analyze the dispersion equation for a ligament solution (55) and carry out the nondimensionalization procedure, the characteristic scales can be selected the same as in the wave solution. Then in a dimensionless form, the dispersion equation (55) is written as:

$$k\_{\ast l}^4 - 2i\varepsilon w\_{\ast} k\_{\ast l}^2 - \varepsilon^2 w\_{\ast}^2 - 3i\varepsilon^2 k\_{\ast l}^3 \omega\_{\ast} - 3\varepsilon^3 k\_{\ast l} w\_{\ast}^2 + \varepsilon k\_{\ast l}^2 - i\varepsilon^2 w\_{\ast} + i\varepsilon^2 w\_{\ast} k\_{\ast l}^3 - \varepsilon^3 k\_{\ast l} w\_{\ast}^2 = 0 \tag{61}$$

Equation (61) has four roots. The analysis shows that the spatial attenuation condition is satisfied for one root only. Due to the cumbersomeness of the expression, the calculated solution is not given here.

Let us construct the dispersion dependences of the components of periodic motion in log–log scales in dimensional variables for liquids with glycerin and water parameters. Figure 4 shows the dependence of the wavelength *λ* on the frequency for the wave component (a) and the analog of the wavelength *δ<sup>i</sup>* on the frequency for the ligament component (b). The letter (*W*) indicates the dependencies for water and the letter (*Gl*) indicates the dependencies for glycerin.

**Figure 4.** Wavelength dependences of the wave solution (**a**) and the ligament solution (**b**) on the frequency in a viscous homogeneous fluid. The curves indicated by the letter (*W*) are constructed for a liquid with water parameters, and by the letter (*Gl*)—for a liquid with glycerin parameters.

Figure 5 shows the dependences of the phase and group velocity on the frequency of periodic motion for the wave component of periodic motion (a) and the ligament associated with the wave component (b). The dependences are constructed for liquids with the parameters of water and glycerin. Figure 6 shows similar dependencies, but on the wavelength. There are several remarkable velocities in a viscous homogeneous liquid. The minimum group velocity and the velocity at which the group and phase velocities are compared. We show that the velocities are compared when the value of the phase velocity is minimal. The extremum condition for the phase velocity will be written as:

$$
\partial\_{\omega}c\_{ph} = \partial\_{\omega} \left(\frac{\omega}{k}\right) = \frac{1}{k^2}(k - \omega \partial\_{\omega}k) = 0\tag{62}
$$

**Figure 5.** Dependences of the phase (dashed lines) and group (solid lines) velocities on the frequency of the wave solution (**a**) and the ligament solution (**b**) on the frequency in a viscous homogeneous fluid. The curves indicated by the letter (*W*) are constructed for a liquid with water parameters, and by the letter (*Gl*) for a liquid with glycerin parameters.

**Figure 6.** Dependences of the phase (dashed lines) and group (solid lines) velocities on the wavelength of the wave solution (**a**) and the ligament solution (**b**) on the frequency in a viscous homogeneous fluid. The curves indicated by the letter (*W*) are constructed for a liquid with water parameters, and by the letter (*Gl*) for a liquid with glycerin parameters.

The ratio (62) leads us to equality:

$$\mathcal{c}\_{ph} = \frac{\omega}{k} = \left(\partial\_{\omega}k\right)^{-1} = \mathcal{c}\_{\mathbb{S}^{\mathbb{T}}}\tag{63}$$

The minimum value of the group velocity for a liquid with water parameters is *cW gr*min = 17.71 cm/s and is achieved at frequency <sup>ω</sup>*<sup>W</sup> gr*min = 40.53 *<sup>s</sup>*−<sup>1</sup> and wavelength *λ<sup>W</sup> gr*min = 4.33 cm. For a liquid with glycerin parameters, the corresponding values are as follows: *cGl gr*min = 29.39 cm/s at the frequency <sup>ω</sup>*<sup>W</sup> gr*min = 21.05 *<sup>s</sup>*−<sup>1</sup> and wavelength *λGl gr*min = 14.88 cm. The minimum phase velocity for a liquid with water parameters is *c<sup>W</sup> ph*min = 23.05 cm/s at the frequency <sup>ω</sup>*<sup>W</sup> ph*min = 84.82 <sup>s</sup>−<sup>1</sup> and wavelength *<sup>λ</sup><sup>W</sup> ph*min = 1.71 cm. For a liquid with glycerin parameters, the corresponding values are as follows: *cGl ph*min = 40.85 cm/s at frequency <sup>ω</sup>*<sup>W</sup> ph*min = 38.96 s−<sup>1</sup> and wavelength *<sup>λ</sup>Gl ph*min = 6.59 cm.

The viscosity of the liquid affects the capillary wave motion. An increase in viscosity leads to an increase in the wavelength at a constant frequency in the region of capillary waves. As a consequence, the values of the phase and group velocity increase. The characteristic values of fluid velocities for the wave component of periodic motion increase with increasing viscosity and shift to the region of lower frequencies (longer wavelengths). Taking into account the viscosity in the model makes it possible to calculate, in addition to the wave component, the ligament component of periodic motion in a liquid.

#### **4. Reduction to Inviscid Fluid**

In an inviscid exponentially stratified fluid, the equations of motion and boundary conditions (12), (13) are reduced ν → 0:

$$z < \zeta : \begin{cases} \rho = \rho\_{00}(r(z) + s(x, z, t)) \\ \rho \partial\_t \mathbf{u} + \rho(\mathbf{u} \cdot \nabla) \mathbf{u} = \rho \mathbf{g} - \nabla P \\ \partial\_t \rho + \text{div}(\rho \mathbf{u}) = 0 \end{cases} \tag{64}$$

$$z < 0 \colon \partial\_{tt} \Delta \psi + N^2 \exp(-z/\Lambda) \partial\_{xx} \psi = 0 \tag{65}$$

$$z = 0 \colon \begin{cases} \partial\_t \zeta + \partial\_x \psi = 0 \\ -\rho \zeta \zeta + \mu \eta h o\_{00} P + \sigma \partial\_{xx} \zeta = 0 \end{cases} \tag{66}$$

Dispersion relations in an inviscid liquid allow us to obtain only the wave component of periodic motion. The dispersion relations for the wave component are obtained using the limit transition ν → 0 from expressions (16) and (30) for high frequency (compared to the buoyancy frequency) and low frequency oscillations, respectively. For high-frequency wave disturbances of the free surface (ω > *N*) we obtain:

$$\begin{cases} -N^2 k\_x^2 + \epsilon^{z/\Lambda} \left(k\_x^2 - k\_z^2\right) \omega^2 = 0\\ \lg k\_x^2 - k\_z \omega^2 + k\_x^4 \gamma = 0 \end{cases} \tag{67}$$

For low frequency waves (ω < *N*):

$$\begin{cases} -N^2 k\_x^2 + \varepsilon^{z/\Lambda} \left(k\_x^2 + k\_z^2\right) \omega^2 = 0\\ -N^2 k\_x^2 - \varepsilon^{z/\Lambda} \left(k\_s^2 - k\_x^2\right) \omega^2 = 0 \end{cases} \tag{68}$$

Thus, the approximation of an ideal fluid enables us to find solutions that describe the wave constituent of periodic motion in the region of gravitational waves quite well. In the region of capillary waves, the wave constituent calculated in the model of an ideal liquid underestimates the value of the wavelength at a given frequency. Note that in the model of an ideal fluid, it is impossible to obtain a solution for the ligament constituent of periodic motion along the free surface of the fluid. This makes the solution incomplete.

#### **5. Discussion**

For the first time the theory of singular perturbations was used to describe the propagation of two-dimensional periodic perturbations over the surface of a viscous exponentially stratified incompressible fluid occupying the lower half-space. Calculations were performed in a single formulation in a wide frequency range, which includes the propagation of capillary, gravitational and internal waves. The system of equations, which includes the equations of continuity, momentum transfer and density profile with depth, replacing the equation of state, is solved together with physically justified kinematic and dynamic conditions on a free surface, taking into account the effect of surface tension. The theory of singular perturbations with respect to the compatibility condition was applied to study

the propagation of infinitesimal harmonic flows with a positive definite frequency and a complex wavenumber, which considers spatial attenuation.

A complete set of solutions of the obtained dispersion relation for infinitesimal periodic perturbations includes two inseparable constituents of periodic flows. The theory of regular perturbations determines the parameters of gravitational-capillary or infra-low-frequency waves with a frequency lower than the buoyancy frequency.

The accompanying fine constituents that are ligaments are calculated by the methods of the theory of singular perturbations. The dispersion equation designed for infinitesimal periodic perturbations describes the relationship between the components of the wave vector in periodic flows in a wide frequency range 10−<sup>4</sup> < ω < 103 s−<sup>1</sup> which includes infra-low frequency, gravitational and capillary waves. The result, represented by lengthfrequency dependence, is convenient for experimental determination. The dependences of phase and group wave velocities on frequency and wavelength have also been obtained.

The dependences of phase and group velocities of regular perturbations (i.e., waves in water and glycerin) on frequency and wavelength have been represented in graphs. Wave properties are significantly modified during the transition through buoyancy frequency. The group velocity of wave propagation tends to zero and the phase one tends to infinity with the approximation of wave and buoyancy frequencies. Viscosity has a significant influence on short waves when their length becomes compared to or less than capillary scale *λ* < *δ γ <sup>g</sup>* = ,*γ*/*g*.

The wave parameters are calculated in an actual homogeneous liquid in the limit of the buoyancy frequency *N* → 0 as well as in a constant density approximation ρ ≡ const, *Tb* = 0 taking into account the compatibility conditions. These calculations correspond to the given results.

Singular solutions describe ligaments (i.e., thin currents accompanying waves in a viscous stratified or homogeneous liquid). To compare them we have shown the graphs of the dependence of the wavelength and scale of the ligaments on the wave frequency in a viscous homogeneous liquid. The sizes of the ligaments differ by several orders of magnitude at low frequencies. The dependences of the phase and group velocities of waves and ligaments on the frequency and wavelength are also given.

The application of the theory of singular perturbations makes it possible to design complete solutions of the dispersion equation and the system of fundamental equations without additional hypotheses and limitations of the type of boundary layer approximation [33,34].

The evolving approach admits the extrapolation to the investigation of three-dimensional periodic perturbations. In this case, propagating surface waves accompany several types of ligaments, as well as internal waves in the thickness of a continuously stratified liquid [68,69].

In general cases, both waves and ligaments continuously interact with each other and generate new groups of flow constituents [75]. The influence of forcing together with the effects of nonlinearity and dissipation causes the evolution of the wave structure, which has different shapes at the phases of growth and attenuation.

In the low frequency range, the forcibly oscillating free surface of a stratified ocean can be a source of internal and inertial waves that transfer energy into the ocean. Internal waves, in turn, interact non-linearly with each other [76], generating new wave groups and ligaments (i.e., accompanying flows that cause the observed modulation of the surface waves) [77]).

#### **6. Conclusions**

The analysis of the linearized reduced version of the fundamental equations system by the methods of the theory of singular perturbations has been carried out with respect to the compatibility condition. It showed that the complete solutions describe waves propagating over the surface of a viscous stratified incompressible fluid, and small-scale constituents that are ligaments accompanying the waves. In extreme cases that are in limit of viscous homogeneous liquid, ideal stratified and ideal homogeneous liquid, the obtained dispersion relations for waves transfer to the widely known ones.

The experimental studies of the fine structure of surface waves in a continuously stratified liquid with high-resolution instruments that allow recording the influence of all constituents of periodic flows are of great interest. We should pay particular attention to the use of high-resolution optical methods for recording variations in the density gradient, and highly sensitive temperature or electrical conductivity sensors in a liquid with salt stratification.

**Author Contributions:** Conceptualization, Y.D.C. and A.A.O.; methodology, Y.D.C. and A.A.O.; software, A.A.O.; validation, Y.D.C. and A.A.O.; formal analysis, Y.D.C.; investigation, Y.D.C.; resources, Y.D.C.; data curation, Y.D.C.; writing—original draft preparation, Y.D.C. and A.A.O.; writing—review and editing, Y.D.C.; visualization, A.A.O.; supervision, Y.D.C.; project administration, Y.D.C.; funding acquisition, Y.D.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by the Russian Science Foundation (project 19-19-00598" Hydrodynamics and energetics of drops and droplet jets: formation, motion, break-up, interaction with the contact surface" https://rscf.ru/en/project/19-19-00598/. (accessed on 12 August 2022).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


*Article*
