Periodic Flows in a Viscous Stratified Fluid in a Homogeneous Gravitational Field

1. Introduction

In natural, laboratory and industrial conditions, the density of a liquid or gas depends on temperature, pressure, concentration of dissolved substances or suspended particles. It is not a constant and changes under the influence of a large number of physical factors. An oscillating source forms waves that propagate over long distances. Historically, it is customary to distinguish acoustic waves, the existence of which is due to the compressibility of the medium, gravitational waves associated with the action of gravity, inertial waves in a globally rotating medium, capillary waves at the interface between media and a large group of hybrid waves, the existence of which is due to the combined action of a number of factors [1, 2].

In the mass forces (gravity, inertia) field the fluid medium is separated. Heavy particles sink, light particles float up - the medium is naturally stratified. Compressibility under the action of hydrostatic pressure is an additional impact on density. The choice of the coordinate system depends on the overall geometry of the problem. Under the Earth conditions, consideration of flows with scales much smaller than its radius is carried out in a Cartesian coordinate system with an axis $z$ pointing vertically upwards. The acceleration of gravity $g$ is directed downward.

The density distribution in the gravity direction $\rho \left(z\right)$ is characterized by the scale $\Lambda ={\left|d\ln \rho \left(z\right)/dz\right|}^{-1}$, frequency $N=\sqrt{g/\Lambda }$ and period $T_b=2\pi /N$ of buoyancy. In the atmosphere and ocean, the average buoyancy period lies in the range of $3<T_b<10$ min [3-5]. In the "instantaneous" density profiles of the atmosphere and ocean, thin, highly gradient interfaces are expressed, separating thick, more homogeneous layers forming a “fine structure” of the medium [3,5].

In practice, several characteristic types of average density distributions have been identified. Further, models of continuous (linear or exponential), two-layer or multi-layer stratification (the last two with a persistent density gap) will be used. In a large group of flows, density variations are much less than its average value.

Observing fluctuations of the free surface and the interface between water and oil in a ship's lighting lamp at sea and mounted on a swing, B. Franklin at the end of the 18^th century noted the need to analyze the influence of fluid density heterogeneity in mathematical research [6]. Initially, the effects of stratification began to be taken into account in calculations of the propagation of internal waves in the atmosphere and ocean, which were carried out by famous English scientists - G.G. Stokes [7], Lord Rayleigh [8], H. Lamb [9] and others.

A systematic study of the influence of stratification on the pattern of flows in the atmosphere and ocean, navigation (the “dead water” effect), noticed in ancient times, began after the publication of the scientific results of F. Nansen’s Polar Expeditions [10]. V. Ekman developing the methodology and planning experiments. In order to conduct laboratory studies of the phenomenon of "dead water" he used a review of the first publications on the theory of internal waves in the treatise [11]. In a series of thorough experiments, V. Ekman determined the conditions for the generation of large waves by a moving model of a ship at a smoothed interface between fresh and salt (sea) water, and determined the influence of the movement mode on the position of the model’s hull and resistance [10]. However, in general, the work on the consideration of the equations of internal waves and the “exotic” phenomenon of “dead water” fell out of scientific circulation for more than half a century and did not affect the development of the general theory of fluid flows.

Among the reasons, at least two have to be noted: the smallness of the density variations compared to its average value, limiting the effect on inertial properties, and the insufficient development of the mathematical apparatus. As G.G. Stokes noted in a fundamental article [12], written several years before a thorough study of wave propagation in homogeneous and layer-by-layer stratified media [13] “As it is quite useless to consider cases of the utmost degree of generality, I shall suppose the fluid to be homogeneous...” However, a few years later, he also emphasized the limitations of the approximation used: “The three equations of which (l) is the type are not the general equations of motion which apply to a heterogeneous fluid when internal friction is taken into account, which are those numbered (10) in my former paper, but are applicable to a homogeneous incompressible fluid, or to a homogeneous elastic fluid subject to small variations of density, such as those which accompany sonorous vibrations” [7].

Accordingly, when studying waves of other types - acoustic [14] or gravitational-capillary at the interface between the atmosphere and the hydrosphere [1, 15], the unperturbed density was assumed to be homogeneous. Here and further, general rotation effects and associated inertial waves [1, 16] will not be considered.

Interest to the mathematical study of the stratification influence began to form in the middle of the last century. In this period precision instruments identified the thin highly gradient structure of the waters of the Baltic Sea [5]. Next, flows induced by diffusion on an inclined wall in a continuously stratified atmosphere were discovered [17]. The development of interest in studying the influence of stratification was facilitated by the papers [18, 19], which showed the important role of diffusion induced flows on topography not only in the atmosphere, where they manifest themselves in the form of mountain and valley winds, but also in the ocean. At the same time, experimental [20] and theoretical studies of internal waves in continuously stratified media [21] began to develop. Numerous expeditions have shown the existence of fine structure and its influence on the dynamics of the atmosphere and ocean in various regions of the Earth.

The number of original articles and reviews describing the influence of stratification on individual phenomena (internal waves, currents, vortices) began to increase rapidly. The propagation of acoustic vibrations in a continuously stratified medium was considered [22]. The influence of viscosity, which was initially taken into account only in terms of exponential attenuation of wave amplitudes [1, 21], began to be analyzed in more detail when describing the propagation of gravitational surface [23-29], internal [30] and acoustic waves [31], taking into account the boundary layers formed simultaneously with the waves.

From the general content of papers and monographs [1,2,11,15,16], it follows that the basis of a rational mathematical description of inhomogeneous fluid flows is a system of fundamental equations - differential analogues of the momentum, energy and matter conservation laws with physically justified initial and boundary conditions. All the equations that were first presented in the first edition of the treatise [1], published in 1944, are quite complex for general analysis. In practice, reduced forms of the general system of equations are usually used, which make it possible to study the properties of individual components of flows - waves, vortices, jets, wakes with the required degree of completeness. In this work, the main attention will be paid to the analysis of periodic flows, the temporal variability of which is proportional to a function of the form $f\propto \exp \left(-i\omega t\right)$.

In the experiment, as at the early stage of development of the analytical theory of waves [11], it was emphasized that the measured physical quantities - parameters of periodic flows, such as the period $T_w$ (frequency $\omega ,$), length $\lambda$, group $c_g$ and phase $c_{ph}$ velocity of the wave are characterized by real numbers. From the very beginning of theoretical study, periodic flows began to be described by complex numbers, introduced to reduce notation and convenience of calculations. Immersion of problems in the algebra of complex numbers leads to an expansion of the dimension of the problem space and the emergence of additional “physically unrealizable” solutions. Accordingly, there is a need to select a part of the solutions corresponding to the initial formulation, with the introduction of criteria explaining the procedure.

The physical interpretation of the solutions depends on the choice of the algorithm for the rules for immersing the problem in the algebra of complex numbers. Traditionally, starting with the works of scientists of the 19^th century, the frequency $\omega$ of a waveform $f\propto \exp \left(ikx-i\omega t\right)$ is chosen as a complex value. Its real part determines the dispersion relation, the functional relationship between frequency $\omega$ and wave vector $k$, and the imaginary part determines the stability condition and the wave attenuation coefficient [1]. An innumerable number of works, including popular monographs, are devoted to the study of the stability of flows and waves [32, 33]. The history of the development of studies of flow stability is traced in detail in [34]. Researchers consider a problem of finding the shape of the liquid surface and the criteria for the development of instability under the action of various destabilizing factors: surface electric charge (Tonks-Frenkel instability) [35, 36], Rayleigh-Taylor and Marangoni thermal convective instabilities [37], etc.

However, from a consideration of the experimental patterns of non-dispersive waves propagation in a medium at rest, it follows that the amplitude and wavelength change with distance from the source, but the frequency of periodic motion remains constant. In this regard, it is natural to keep the frequency, a measure of wave energy, as a positive definite real quantity in calculations, and take the wave number to be complex [38]. Substituting expansions of this type into a linearized system of fundamental equations, the solution of which is found by methods of singular perturbation theory [39], allows for a new classification of the structural components of periodic flows based on the properties of complete solutions.

That part of the solutions of the fundamental equations system, that includes regularly perturbed functions, characterizes waves slowly decaying in the direction of propagation in weakly dissipative media. Singularly perturbed components of the solution describe ligaments—thin flows that determine the structure of the medium in both linear and weakly nonlinear approximations [40,41,42].

In the hydrosphere and atmosphere, there are types of waves that differ significantly in frequency (in particular, acoustic and internal waves in the thickness of stratified liquid [1, 11,21,22,31]) or in the distribution of displacement amplitudes in depth (surface and internal waves [1 ,11,21]. That makes it possible to study their properties within the framework of individual specialized equations - acoustics [1,22,31], internal [21] surface gravity or capillary waves [11,43].

The patterns of propagation of a set of two-dimensional periodic disturbances - waves and ligaments - in an incompressible fluid, when the reduced continuity equation allows us to introduce a stream function convenient for analysis, are considered in the thickness [38] and on the surface of a viscous stratified fluid [44]. This paper is the first to consider the problem of propagating a complete set of two-dimensional infinitesimal periodic disturbances in a continuously stratified compressible fluid.

2. System of fundamental equations of periodic flows in the atmosphere and ocean.

2.1. The complete system of equations determining the flow of the liquid.

Periodic wave processes occurring in a viscous liquid are considered. The liquids existing in nature are heterogeneous. The inhomogeneous distribution of density $\rho$ is determined by the equation of state

$$\rho =\rho \left(P,S,s_n,T\right).$$

(1)

Here the symbol $P$ denotes pressure, $S$ stands for entropy, $s_n$ denotes salinity of the n-th impurity and $T$ stands for temperature

Far from the conditions of phase transitions, the values of the temperature gradient and the impurity content are limited and it is permissible to use a linearized equation of state:

$$\begin{array}{c}\rho ={\rho }_0\left(1-{\alpha }_T\left(T-T_0\right)+{\alpha }_P\left(P-P_0\right)+{\displaystyle \sum _n{\alpha }_{s_n}\left(s_n-s_{n0}\right)}\right)\\ {\alpha }_T=-\frac{1}{\rho }\frac{\partial \rho }{\partial T},{\alpha }_P=\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial P}\right)}_S,{\alpha }_{s_n}=\frac{1}{\rho }\frac{\partial \rho }{\partial s_n}\end{array}.$$

(2)

Here ${\alpha }_T$ denotes the coefficient of thermal expansion of the liquid, ${\alpha }_P$ stands for the coefficient of adiabatic compressibility of the liquid, ${\alpha }_{s_n}$ denotes the coefficient of contraction of the n-th impurity and $T_0,P_0,s_{n0}$ stand for reference level of temperature, pressure and salinity respectively.

Fundamental system of equations in addition to the equation of state consists of equations for describing the transfer of matter, concentration of impurity, temperature, and momentum. Taking into account the neglect of the thermophoresis Ludwig–Soret effect [45 – 46] and Dufour effect [47], the system of equations is written as follows [38, 44]:

$$G=G\left(P,S,s_n,T\right),\rho =\rho \left(P,S,s_n,T\right)$$

(3)

$${\partial }_t\rho +{\nabla }_j\left(p^j\right)=Q_{\rho },$$

(4)

$${\partial }_t\left(p^i\right)+{\nabla }_j{\Pi }^{ij}=\rho g^i+2{\epsilon }^{ijk}{p}_j{\mathsf{\Omega }}_k+Q^i,$$

(5)

$${\partial }_t\left(\rho T\right)+{\nabla }_j\left(p^{j}T\right)=\Delta \left({\kappa }_T\rho T\right)+Q_T,$$

(6)

$${\partial }_t\left(\rho s_n\right)+{\nabla }_j\left(p^j{s}_n\right)=\Delta \left({\kappa }_s\rho s_n\right)+Q_{s_n}.$$

(7)

Here $G$ is Gibbs potential, $Q_{\rho },Q^i,Q_T,Q_{s_n}$ stand for the source of mass, momentum, temperature and salinity concentration, respectively, $p$ denotes momentum, ${\Pi }^{ij}=\rho u^i{u}^j+P{\delta }^{ij}-{\sigma }^{ij}$ stands for the momentum flux density tensor, $u^i$ is components of the fluid velocity $u=p/\rho$, ${\delta }^{ij}$ is the Kronecker delta, ${\sigma }^{ij}=\mu \left(\frac{\partial u^i}{\partial x^j}+\frac{\partial u^j}{\partial x^i}-\frac{2}{3}{\delta }^{ij}\frac{\partial u^k}{\partial x^k}\right)+\zeta {\delta }^{ij}\frac{\partial u^k}{\partial x^k}$ denotes the viscous stress tensor, $\mu ,\zeta$ is dynamic and bulk viscosity respectively, $g$ is the gravity acceleration, $\epsilon$ is the internal energy, $\mathsf{\Omega }$ is the global rotation angular velocity and ${\kappa }_T,{\kappa }_{s_n}$ stand for thermal and mass diffusivity respectively.

Equations (1), (3) – (7) form a fundamental system of equations that determine fluid flow. The complete solution of the system of equations (3) – (7) determines all components of flow in liquids - waves: acoustic, gravitational (internal and surface), capillary, hybrid and ligaments - accompanying components that determine the fine structure of the flow. Usually researchers ignore the fine structure, limiting themselves to a partial solution of a system of equations. In this work, we construct a theory that takes into account all flow components.

To complete the formulation, it is necessary to add initial and boundary conditions to the problem. The initial conditions depend on the shape and type of the oscillation source. Often, when studying the properties of periodic flows, instead of initial conditions, they specify the type of solution and look for steady-state solutions of a given type. On the surface of a solid impermeable body $\Sigma$ no-slip, no-flux boundary and initial conditions is written as follows:

$${\left.u\right|}_{\Sigma }=0,{\left.u\right|}_{t\le 0}=0,{\left.P\right|}_{t\le 0}=P_0,{\left.s_n\right|}_{t\le 0}=s_{n0},{\left.T\right|}_{t\le 0}=T_0,$$

(8)

If the distance to the boundaries greatly exceeds the characteristic dimensions of the observed phenomena, then a model of an unbounded medium is often used. In this case, the boundary conditions are transformed into the conditions of physical implementation - attenuation with removal:

$${\left.u\right|}_{r\to \infty }\to 0,$$

(9)

If the model under consideration contains a free surface or interface between layers of immiscible liquids, then it is necessary to add standard hydrodynamic boundary conditions: kinematic and dynamic boundary conditions. The kinematic boundary condition is written for both contacting layers (or for one medium in the case of a free surface): the substantial derivative of the function $F$ defining the shape of the free surface is equal to zero at the boundary:

$$\frac{DF}{Dt}\equiv \frac{\partial F}{\partial t}+\left(u\cdot \nabla \right)F=0,$$

(10)

Dynamic boundary conditions are determined by the balance of forces at the interface (free surface of the liquid):

$$n_{\left(1\right)}^k{\sigma }_{\left(1\right)}^{ik}+n_{\left(2\right)}^k{\sigma }_{\left(2\right)}^{ik}=0,$$

(11)

Here $n$ is the unit normal vector, and the subscript (1), (2) refers to the two contacting media. If the model takes into account the effects of surface tension, then on the right side of (11) it is necessary to take into account Laplace forces as well.

2.2. The reduced system of equations.

The fundamental system of equations is complete and allows one to determine the patterns of changes in basic physical quantities during the propagation of periodic disturbances in continuous media. Since the complete system of equations is of high order and very complex to analyze, it is simplified to study the properties of individual processes. An extremely simplified model in which it is possible to track the dynamics and evolution of the structure of periodic flows takes into account the uneven distribution of density, without indicating the physical nature of the formation of heterogeneity.

The system of equations (3) – (7) is noticeably reduced in the constant temperature model in the absence of impurities in the weakly compressible fluid. The consideration is carried out in a Cartesian coordinate system $Oxyz$ in which the $Oz$ axis is directed against the direction of the gravity acceleration $g$. The $Oxy$ plane determines the position of the reference level. In a weakly compressible viscous fluid, bulk viscosity takes on a zero value. In the absence of mass sources $Q_{\rho }=0$ and under the assumptions made, the reduced system of equations will take the form:

$$\rho \left({\partial }_{t}u+\left(u\cdot \nabla \right)u\right)=\rho \nu \Delta u-\nabla P+\rho g$$

(12)

$${\partial }_t\rho +u\cdot \nabla \rho +\rho \mathrm{div}u=0$$

(13)

$$\rho ={\rho }_0\left(z\right)\left(1+\tilde{\rho }\left(x,y,z,t\right)\right)$$

(14)

The initial stratification ${\rho }_0\left(z\right)$ when describing models is often defined as a linear ${\rho }_0\left(z\right)={\rho }_{00}\left(1-z/\Lambda \right)$ or exponential ${\rho }_0\left(z\right)={\rho }_{00}\exp \left(-z/\Lambda \right)$f unction. The symbol ${\rho }_{00}$ indicates the density value at the reference level $z=0$, and the symbol $\Lambda ={\left|d\ln \rho /dz\right|}^{-1}$characterizes the stratification scale. A stably stratified liquid is characterized by the limiting frequency of its own mechanical vibrations [8] – buoyancy frequency, the square of which is given by:

$$N^2=-\frac{g}{\rho }\frac{d\rho }{dz}$$

(15)

Equation of state (2) under the assumptions made is simplified:

$$\rho ={\rho }_0\left(z\right)\left(1-{\alpha }_P\left(P-P_0\right)\right)$$

(16)

Fluid pressure is represented as the sum of reference level pressure $P$, hydrostatic pressure and perturbation pressure $\tilde{P}$:

$$P=P_0+{\displaystyle {\int }_z^0\rho \left(x,y,\xi ,t\right)gd\xi +\tilde{P}\left(x,y,z,t\right)}$$

(17)

Taking into account the equation of state (16), the definition of the velocity of sound $c^2={\left(\partial P/\partial \rho \right)}_S$ and the definition of pressure (17), relation (15) for the buoyancy frequency takes the form:

$$N^2=\frac{g^2}{c^2}\left(\frac{c_p}{c_V}-1\right),$$

(18)

Here $c_P,c_V$ is the heat capacity at constant pressure and at constant volume, respectively.

The resulting system of equations, despite significant simplifications, qualitatively completely describes periodic flows in viscous inhomogeneous continuous media. The boundary and initial conditions will not change.

In the model under consideration there are intrinsic parameters. These parameters determine the characteristic scales of the flow components and the characteristic times of their observation. A set of kinetic coefficients allows you to form your own parameters. Intrinsic parameters for liquids with the parameters of water and air are presented in Table 1 and Table 2, respectively.

The natural parameters presented in the table have to be supplemented with temporal and spatial scales that do not depend on the level of fluid stratification. Taking into account compressibility, a time scale ${\tau }_c^{\nu }=\nu /c^2$ is added. It takes values for water ${\tau }_c^{\nu }\simeq 4\cdot {10}^{-13},s$ and for air ${\tau }_c^{\nu }\simeq {10}^{-10},s$. Spatial scale ${\delta }_c^{\nu }=\nu /c$ is added. It takes values for water ${\delta }_c^{\nu }\simeq 7\cdot {10}^{-10},m$ and for air ${\delta }_c^{\nu }\simeq 5\cdot {10}^{-8},m$. In viscous liquids (homogeneous and heterogeneous), a capillary-viscous time scale appears ${\tau }_{\nu g}^{\gamma }=\gamma /\nu g$. The symbol $\gamma =\sigma /{\rho }_{00}$ denotes the surface tension coefficient of the liquid $\sigma$normalized to the equilibrium density value ${\rho }_{00}$. For water, the capillary-viscous time scale takes on values ${\tau }_{\nu g}^{\gamma }\simeq 7s$, and for air ${\tau }_{\nu g}^{\gamma }\simeq 400s$. The spatial scale in viscous liquids ${\delta }_g^{\nu }=\sqrt[3]{{\nu }^2/g}$ has the value ${\delta }_g^{\nu }\simeq 5\cdot {10}^{-5}\mathrm{m}$ in water and ${\delta }_g^{\nu }\simeq 3\cdot {10}^{-4}\mathrm{m}$ in air. Both in the model of a viscous and in the model of an inviscid liquid one of the proper parameters is capillary length ${\delta }_g^{\gamma }=\sqrt{\gamma /g}$. For water capillary length takes the value ${\delta }_g^{\gamma }\simeq 3\cdot {10}^{-3}m$and for air it is ${\delta }_g^{\gamma }\simeq 8\cdot {10}^{-2}m$.

Small disturbances of physical quantities (pressure, density, velocity) are often occurs in nature. Let us solve the problem by the decomposition method for a small parameter that plays the role of the amplitude of periodic movements.

3. Periodic flows in the thickness of a uniformly stratified liquid.

3.1. Linearization of the equations system.

Perturbations of the target values (velocity, density and pressure) will be considered small. To obtain dispersion relations, we linearize the system of equations (12) – (14), (16). If we assume that the fluid is exponentially stratified, then in a linear approximation in terms of the amplitude of periodic motion, the reduced system of fundamental equations is written as follows:

$$\left\{\begin{array}{l}{\partial }_t\tilde{\rho }-\frac{w}{\Lambda }+{\partial }_{x}u+{\partial }_{y}v+{\partial }_{z}w=0\\ {\partial }_{t}u-\nu \Delta u+\frac{1}{{\rho }_{00}}{\partial }_x\tilde{P}=0\\ {\partial }_{t}v-\nu \Delta v+\frac{1}{{\rho }_{00}}{\partial }_y\tilde{P}=0\\ {\partial }_{t}w-\nu \Delta w+\frac{1}{{\rho }_{00}}{\partial }_z\tilde{P}+g\tilde{\rho }=0\\ \frac{1}{{\rho }_{00}{c}^2}{\partial }_t\tilde{P}-\frac{wg}{c^2}+{\partial }_{x}u+{\partial }_{y}v+{\partial }_{z}w=0\end{array}\right.$$

(19)

Here $u,v,w$ are the components of the velocity field $u=\left(u,v,w\right)$. We will look for the solution of the equations system (19) in the form of periodic flows $\propto \exp \left(i\omega t\right)$:

$$\left(\begin{array}{l}u\\ v\\ w\\ \tilde{P}\\ \tilde{\rho }\end{array}\right)=\left(\begin{array}{l}{U}_m\\ V_m\\ W_m\\ P_m\\ {{\rm P}}_m\end{array}\right)\exp \left(ikr-i\omega t\right)=\left(\begin{array}{l}{U}_m\\ V_m\\ W_m\\ P_m\\ {{\rm P}}_m\end{array}\right)\exp \left(i{k}_{x}x+i{k}_{y}y+i{k}_{z}z-i\omega t\right)$$

(20)

Here $U_m,V_m,W_m,P_m,{{\rm P}}_m$ are the amplitudes of the corresponding quantities, $k$ is the wave vector, the components of which have the right to be complex values $k_x,k_y,k_z$, and the frequency of periodic motion $\omega$ is considered positive definite.

3.2. Dispersion relation. Classification of flow components.

By substituting the type of solution (20) into the system of equations (19), we obtain a system of algebraic equations. The compatibility condition of the algebraic equations system determines the dispersion relations between the components of the wave vector and the frequency of periodic motion:

$$\begin{array}{c}{D}_{\nu }\left(k\right)\left({\omega }^2{D}_{\nu }^2\left(k\right)-\omega N^2{D}_{\nu }\left(k\right)+c^2{k}_{\perp }^2{N}_c^2-c^2\omega k^2{D}_{\nu }\left(k\right)\right)=0,\\ D_{\nu }\left(k\right)=\omega +i\nu k^2,k^2=k_x^2+k_y^2+k_z^2,k_{\perp }^2=k_x^2+k_y^2,N^2=\frac{g}{\Lambda },N_c^2=N^2-\frac{g^2}{c^2}\end{array}$$

(21)

Dispersion relation (21) coincides with the relation obtained earlier [48] in which the limiting transition to a non-rotating weakly compressible fluid (second viscosity $\zeta \to 0$) was made. It is convenient to find and analyze the regular and singular components of the solution to the dispersion relation (21) in dimensionless variables if one chooses the own scales of the problem as non-dimensional parameters. Intrinsic scales characterize the spatial and temporal dimensions of the observed phenomena (see Table 1, 2). We will choose the inverse buoyancy frequency ${\tau }_b=N^{-1}$ as the time scale, and the viscous wave scale ${\delta }_N^{g\nu }={\left(g\nu \right)}^{1/3}{N}^{-1}$ as the spatial scale. With the selected non-dimensional parameters, the dispersion relation (21) is written as follows:

$$\begin{array}{c}\left(i{k}_{*}^2\epsilon +{\omega }_{*}\right)\left(k_{\perp *}^2\left(\frac{\epsilon }{\eta }-\frac{1}{{\epsilon }^2}\right)+{\omega }_{*}^2{\left(i{k}_{*}^2\epsilon +{\omega }_{*}\right)}^2-{\omega }_{*}\left(i{k}_{*}^2\epsilon +{\omega }_{*}\right)-k_{*}^2{\omega }_{*}\frac{\epsilon }{\eta }\left(i{k}_{*}^2\epsilon +{\omega }_{*}\right)\right)=0,\\ \epsilon =\frac{{\delta }_g^{\nu }}{{\delta }_N^{g\nu }}=\frac{\sqrt{\nu /N}}{{\left(g\nu \right)}^{1/3}{N}^{-1}}=\frac{N{\nu }^{1/3}}{g^{2/3}},\eta =\frac{{\tau }_c^{\nu }}{{\tau }_b}=\frac{N\nu }{c^2}\end{array}$$

(22)

The ratio of the natural parameters of the medium - viscous ${\delta }_g^{\nu }$and viscous wave scales ${\delta }_N^{g\nu }$and the time scales ratio $\eta$ characterize the small parameters of the problem. The dimensionless components of the wave vector and the dimensionless frequency are indicated by the subscript «*». Since at the highest degree of equation (22) there is a small parameter, the equation is singularly perturbed with respect to $k_{*z}$. Consequently, the solution of the form $k_{*z}=k_{*z}\left(k_{*x},k_{*y},{\omega }_{*}\right)$ contains regular and singular components. Solutions of equation (22) are written as follows:

$$k_{*z}=\pm \sqrt{-k_{*\perp }^2+\frac{i{\omega }_{*}}{\epsilon }}$$

(23)

$$k_{*z}=\pm \sqrt{\frac{-\epsilon {\omega }_{*}\left({\omega }_{*}+2k_{*\perp }^2\epsilon \left(i+\eta {\omega }_{*}\right)-i\eta \left(2{\omega }_{*}^2-1\right)\right)-\sqrt{-{\epsilon }^2{\omega }_{*}^2\left(\eta +i{\omega }_{*}\right)+4{\omega }_{*}{k}_{*\perp }^2\left({\epsilon }^3-\eta \right)\left(i+\eta {\omega }_{*}\right)}}{2{\epsilon }^2{\omega }_{*}\left(i+\eta {\omega }_{*}\right)}}$$

(24)

$$k_{*z}=\pm \sqrt{\frac{-\epsilon {\omega }_{*}\left({\omega }_{*}+2k_{*\perp }^2\epsilon \left(i+\eta {\omega }_{*}\right)-i\eta \left(2{\omega }_{*}^2-1\right)\right)+\sqrt{-{\epsilon }^2{\omega }_{*}^2\left(\eta +i{\omega }_{*}\right)+4{\omega }_{*}{k}_{*\perp }^2\left({\epsilon }^3-\eta \right)\left(i+\eta {\omega }_{*}\right)}}{2{\epsilon }^2{\omega }_{*}\left(i+\eta {\omega }_{*}\right)}}$$

(25)

In dimensional form, roots (23) – (25) will be written as follows:

$$k_z=\pm \sqrt{-k_{\perp }^2+\frac{i\omega }{\nu }}$$

(26)

$$k_{*z}=\pm \sqrt{\frac{-\left(i\nu N^2+2\nu \omega \left(\nu k_{\perp }^2-i\omega \right)+c^2\left(2i\nu k_{\perp }^2+\omega \right)\right)-\sqrt{-{\nu }^2{N}^4+2\nu c^2\left(2\nu k_{\perp }^2{N}_c^2-i\omega N^2\right)+c^4\left({\omega }^2+\frac{4i\nu N_c^2{k}_{\perp }^2}{\omega }\right)}}{2\nu \left(i{c}^2+\nu \omega \right)}}$$

(27)

$$k_{*z}=\pm \sqrt{\frac{-\left(i\nu N^2+2\nu \omega \left(\nu k_{\perp }^2-i\omega \right)+c^2\left(2i\nu k_{\perp }^2+\omega \right)\right)+\sqrt{-{\nu }^2{N}^4+2\nu c^2\left(2\nu k_{\perp }^2{N}_c^2-i\omega N^2\right)+c^4\left({\omega }^2+\frac{4i\nu N_c^2{k}_{\perp }^2}{\omega }\right)}}{2\nu \left(i{c}^2+\nu \omega \right)}}$$

(28)

The choice of sign in solutions (23) – (25) or (26) – (28) is determined by the boundary conditions for the decay of periodic motion with distance from the source of disturbances. Solutions (25) and (28) describe the regular component and determine the wave motion, solutions (23) - (24) and (26) - (27) determine the singular component of the solution and determine two types of ligaments.

Let us consider the behavior of dispersion relations in limiting cases.

4. High-frequency acoustic waves

Let us consider the limit of high-frequency oscillations corresponding to acoustic oscillations if their oscillation frequency significantly exceeds the buoyancy frequency of the medium $\omega \gg N$ [22, 48]. In this approximation, dispersion equation (21) will be rewritten as:

$$D_{\nu }\left(k\right)\left(D_{\nu }\left(k\right)\omega \left(D_{\nu }\left(k\right)\omega -c^2{k}^2\right)-g^2{k}_{\perp }^2\right)=0,$$

(29)

The solution to dispersion relation (29) is written as:

$$\begin{array}{l}{k}_z=\pm \sqrt{-k_{\perp }^2+\frac{i\omega }{\nu }};\\ k_z=\pm \sqrt{-k_{\perp }^2-\frac{c^2\omega -2i\nu {\omega }^2+\sqrt{c^4{\omega }^2-\frac{4g^2\nu k_{\perp }^2}{\omega }\left(i{c}^2+\nu \omega \right)}}{2\nu \left(i{c}^2+\nu \omega \right)}};\\ k_z=\pm \sqrt{-k_{\perp }^2-\frac{c^2\omega -2i\nu {\omega }^2-\sqrt{c^4{\omega }^2-\frac{4g^2\nu k_{\perp }^2}{\omega }\left(i{c}^2+\nu \omega \right)}}{2\nu \left(i{c}^2+\nu \omega \right)}};\end{array}$$

(30)

The roots (30) describe the wave motion and two attached ligaments. The sign in solution (30) is chosen based on the need for attenuation of periodic motion $\mathrm{Im}\left(k_z\right)>0$ when moving in the positive direction of the axis $Oz$. For oppositely directed motion, the solutions are symmetrical.

When moving to a 2D formulation (if we consider the movement to be independent of the horizontal coordinate $y$), one of the ligaments degenerates and the solution contains one wave and one ligament component:

$$\begin{array}{l}{k}_z=\pm \sqrt{-k_x^2-\frac{c^2\omega -2i\nu {\omega }^2+\sqrt{c^4{\omega }^2-\frac{4g^2\nu k_x^2}{\omega }\left(i{c}^2+\nu \omega \right)}}{2\nu \left(i{c}^2+\nu \omega \right)}};\\ k_z=\pm \sqrt{-k_x^2-\frac{c^2\omega -2i\nu {\omega }^2-\sqrt{c^4{\omega }^2-\frac{4g^2\nu k_x^2}{\omega }\left(i{c}^2+\nu \omega \right)}}{2\nu \left(i{c}^2+\nu \omega \right)}};\end{array}$$

(31)

In the limit of an inviscid fluid, dispersion relation (29) is simplified even further and written in form:

$${\omega }^2\left({\omega }^2-k^2{c}^2\right)-g^2{k}_{\perp }^2=0,$$

(32)

The ligament components of the solution to relation (32) degenerate and only the wave component remains:

$$k_z=\pm \sqrt{-k_{\perp }^2-\frac{g^2{k}_{\perp }^2}{c^2{\omega }^2}+\frac{{\omega }^2}{c^2}}$$

(33)

5. Low-frequency gravity waves

In the limit of low-frequency oscillations $\omega \ll N$, the dispersion relation (21) takes the form:

$$D_{\nu }\left(k\right)\left(c^2\omega i{k}^4\nu -c^2{N}^2{k}_{\perp }^2+c^2{k}^2{\omega }^2+N^2\omega D_{\nu }\left(k\right)+g^2{k}_{\perp }^2\right)=0,$$

(34)

Relation (34) also contains a solution in the form of a wave disturbance and two attached ligaments:

$$\begin{array}{l}{k}_z=\pm \sqrt{-k_{\perp }^2+\frac{i\omega }{\nu }};\\ k_z=\pm \sqrt{-k_{\perp }^2-\frac{i{c}^2{\omega }^2-N^2\nu \omega +\sqrt{N^4{\nu }^2{\omega }^2+4i{c}^2{k}_{\perp }^2\omega \nu \left(g^2-c^2{N}^2\right)-c^4{\omega }^4}}{2c^2\nu \omega }};\\ k_z=\pm \sqrt{-k_{\perp }^2-\frac{i{c}^2{\omega }^2-N^2\nu \omega -\sqrt{N^4{\nu }^2{\omega }^2+4i{c}^2{k}_{\perp }^2\omega \nu \left(g^2-c^2{N}^2\right)-c^4{\omega }^4}}{2c^2\nu \omega }};\end{array}$$

(35)

When transitioning to a flat formulation, one of the ligaments degenerates:

Relation (34) also contains a solution in the form of a wave disturbance and two attached ligaments:

$$\begin{array}{l}{k}_z=\pm \sqrt{-k_x^2-\frac{i{c}^2{\omega }^2-N^2\nu \omega +\sqrt{N^4{\nu }^2{\omega }^2+4i{c}^2{k}_x^2\omega \nu \left(g^2-c^2{N}^2\right)-c^4{\omega }^4}}{2c^2\nu \omega }};\\ k_z=\pm \sqrt{-k_x^2-\frac{i{c}^2{\omega }^2-N^2\nu \omega -\sqrt{N^4{\nu }^2{\omega }^2+4i{c}^2{k}_x^2\omega \nu \left(g^2-c^2{N}^2\right)-c^4{\omega }^4}}{2c^2\nu \omega }};\end{array}$$

(36)

In an ideal liquid, dispersion relation (34) is simplified:

$$c^2{k}^2{\omega }^2-c^2{k}_{\perp }^2{N}^2+N^2{\omega }^2+g^2{k}_{\perp }^2=0,$$

(37)

Solution (37), which is represented only by a wave component, in an ideal liquid transforms into a well-known expression that does not include the wavelength ${\omega }^2={\sin }^2\mathsf{\theta }\cdot N^2$, which describes the geometry of the wave packet in the shape of a “St. Andrew’s cross” ($\mathsf{\theta }$ is the angle of inclination of the wave vector to the horizontal) [20, 21 ]. Ligaments in the ideal fluid model degenerate:

$$k_z=\pm \sqrt{-k_{\perp }^2+k_{\perp }^2\frac{N^2{c}^2-g^2}{c^2{\omega }^2}-\frac{N^2}{c^2}};$$

(38)

The limiting cases discussed in paragraphs 4 and 5 show that ligaments are observed in the entire frequency range from infra-low-frequency mechanical vibrations to high-frequency sound vibrations. The fine structure of the flow accompanies wave motion and requires attention when analyzing phenomena.

6. Periodic flows in two-layer system of stratified liquids

In a two-layer system, a stratified weakly compressible ocean and a stratified compressible atmosphere, it is necessary to write down the boundary conditions at the interface. In a two-layer system, the pressure in both media is written in the form of the sum of hydrostatic pressure and perturbation pressure and in a 2D formulation (if we consider the movement independent of the horizontal coordinate $y$) it is written as follows:

$$P^{o,a}={\displaystyle {\int }_z^{\zeta }{\rho }^{o,a}\left(x,\xi ,t\right)gd\xi +{\tilde{P}}^{o,a}\left(x,z,t\right)}$$

(39)

Here and further, the superscripts “o” and “a” denote quantities related to the ocean (the lower denser liquid) and the atmosphere (the upper less dense liquid), respectively. The symbol $\zeta =\zeta \left(x,t\right)$denotes the function that determines the deviation of the interface between media from the equilibrium position $z=0$. The system of equations of motion, taking into account expressions (39) is written as follows:

$$z<\zeta :{\partial }_t{u}^o-{\nu }^o\Delta u^o+\frac{1}{{\rho }_{00}^o}\nabla P^o-{\rho }^{o}g=0$$

(40)

$${\partial }_t{\rho }^o+u^o\cdot \nabla {\rho }^o+{\rho }^o\mathrm{div}{u}^o=0$$

(41)

$${\rho }^o={\rho }_0^o\left(z\right)\left(1-{\alpha }_P^o\left(P^o-P_0^o\right)\right)$$

(42)

$$z>\zeta :{\partial }_t{u}^a-{\nu }^a\Delta u^a+\frac{1}{{\rho }_{00}^a}\nabla P^a-{\rho }^{a}g=0$$

(43)

$${\partial }_t{\rho }^a+u^a\cdot \nabla {\rho }^a+{\rho }^a\mathrm{div}{u}^a=0$$

(44)

$${\rho }^a={\rho }_0^a\left(z\right)\left(1-{\alpha }_P^a\left(P^a-P_0^a\right)\right)$$

(45)

The system of equations (40) – (45) is supplemented with boundary conditions at the interface: $z=\xi$

$$z=\zeta :{\partial }_t\zeta +u^o{\partial }_x\zeta =w^o$$

(46)

$${\partial }_t\zeta +u^a{\partial }_x\zeta =w^a$$

(47)

$$P^o-2{\rho }^o{\nu }^{o}n\cdot \left(\left(n\cdot \nabla \right)u^o\right)=P^a-2{\rho }^a{\nu }^{a}n\cdot \left(\left(n\cdot \nabla \right)u^a\right)-\sigma \mathrm{div}n$$

(48)

$$u^o\cdot \tau =u^a\cdot \tau$$

(49)

$${\rho }^o{\nu }^o\left(\tau \cdot \left(\left(n\cdot \nabla \right)u^o\right)+n\cdot \left(\left(\tau \cdot \nabla \right)u^o\right)\right)={\rho }^a{\nu }^a\left(\tau \cdot \left(\left(n\cdot \nabla \right)u^a\right)+n\cdot \left(\left(\tau \cdot \nabla \right)u^a\right)\right)$$

(50)

$$n=\frac{\nabla \left(z-\zeta \right)}{\left|\nabla \left(z-\zeta \right)\right|}=\left(\frac{-{\partial }_x\zeta }{\sqrt{1+{\left({\partial }_x\zeta \right)}^2}},\frac{1}{\sqrt{1+{\left({\partial }_x\zeta \right)}^2}}\right),\tau =\left(\frac{1}{\sqrt{1+{\left({\partial }_x\zeta \right)}^2}},\frac{{\partial }_x\zeta }{\sqrt{1+{\left({\partial }_x\zeta \right)}^2}}\right)$$

Here $\sigma$ is the coefficient of surface tension at the interface between contacting media, and $n,\tau$ are the normal and tangent vectors to the interface respectivetly. After carrying out the linearization procedure and transferring the boundary conditions to the equilibrium surface $z=0$ [49] the mathematical formulation in a linear approximation takes the form:

$$\begin{array}{c}z<0:{\displaystyle {\int }_z^{\zeta }{\mathrm{e}}^{-\frac{\xi }{{\Lambda }^o}}g{\partial }_x{\tilde{\rho }}^o\left(x,\xi ,t\right)+{\mathrm{e}}^{-\frac{\zeta }{{\Lambda }^o}}g{\partial }_x\zeta +{\partial }_t{u}^o-{\nu }^o\Delta u^o+{\partial }_x{\tilde{P}}^o}=0\\ {\mathrm{e}}^{-\frac{z}{{\Lambda }^o}}{\partial }_t{w}^o-{\nu }^o{\mathrm{e}}^{-\frac{z}{{\Lambda }^o}}\Delta w^o+\frac{{\partial }_z{\tilde{P}}^o}{{\rho }_{00}^o}=0\\ {\partial }_t{\tilde{\rho }}^o-\frac{w^o}{{\Lambda }^o}+{\partial }_x{u}^o+{\partial }_z{w}^o=0\\ \frac{1}{{\rho }_{00}^o{c}^{o2}}{\partial }_t{\tilde{P}}^o-\frac{w^{o}g}{c^{o2}}+{\partial }_x{u}^o+{\partial }_z{w}^o=0\end{array}$$

(51)

$$\begin{array}{c}z>0:{\displaystyle {\int }_z^{\zeta }{\mathrm{e}}^{-\frac{\xi }{{\Lambda }^a}}g{\partial }_x{\tilde{\rho }}^a\left(x,\xi ,t\right)+{\mathrm{e}}^{-\frac{\zeta }{{\Lambda }^a}}g{\partial }_x\zeta +{\partial }_t{u}^a-{\nu }^a\Delta u^a+{\partial }_x{\tilde{P}}^a}=0\\ {\mathrm{e}}^{-\frac{z}{{\Lambda }^a}}{\partial }_t{w}^a-{\nu }^a{\mathrm{e}}^{-\frac{z}{{\Lambda }^a}}\Delta w^a+\frac{{\partial }_z{\tilde{P}}^a}{{\rho }_{00}^a}=0\\ {\partial }_t{\tilde{\rho }}^a-\frac{w^a}{{\Lambda }^a}+{\partial }_x{u}^a+{\partial }_z{w}^a=0\\ \frac{1}{{\rho }_{00}^a{c}^{a2}}{\partial }_t{\tilde{P}}^a-\frac{w^{a}g}{c^{a2}}+{\partial }_x{u}^a+{\partial }_z{w}^a=0\end{array}$$

(52)

$$\begin{array}{c}z=0:{\partial }_t\zeta -w^o=0,{\partial }_t\zeta -w^a=0,u^o-u^a=0,\\ {\tilde{P}}^o-{\tilde{P}}^a+2{\rho }_{00}^a{\nu }^a{\partial }_z{w}^a-2{\rho }_{00}^o{\nu }^o{\partial }_z{w}^o+\sigma {\partial }_{xx}\zeta =0,\\ {\rho }^o{\nu }^o\left({\partial }_z{u}^o+{\partial }_x{w}^o\right)-{\rho }^a{\nu }^a\left({\partial }_z{u}^a+{\partial }_x{w}^a\right)=0\end{array}$$

(53)

We will look for a solution to the system of equations (51) – (53) in the form of periodic flows $\propto \exp \left(i\omega t\right)$:

$$\left(\begin{array}{l}{u}^{o,a}\\ w^{o,a}\\ {\tilde{P}}^{o,a}\\ {\tilde{\rho }}^{o,a}\\ \zeta \end{array}\right)=\left(\begin{array}{l}{U}_m^{o,a}\\ W_m^{o,a}\\ P_m^{o,a}\\ {{\rm P}}_m^{o,a}\\ A_m\end{array}\right)\exp \left(i{k}^{o,a}r-i\omega t\right)=\left(\begin{array}{c}{U}_m^{o,a}\exp \left(i{k}_z^{o,a}z\right)\\ W_m^{o,a}\exp \left(i{k}_z^{o,a}z\right)\\ P_m^{o,a}\exp \left(i{k}_z^{o,a}z\right)\\ {{\rm P}}_m^{o,a}\exp \left(i{k}_z^{o,a}z\right)\\ A_m\end{array}\right)\exp \left(i{k}_{x}x-i\omega t\right)$$

(54)

Substituting the type of solution (54) into the main equations (51) – (52) leads to a system of algebraic equations connecting the components of wave vectors $k_x,k_z^{o,a}$ and the frequency of periodic disturbances $\omega$:

$$\left(\begin{array}{cccc}-g{k}_x^2+\frac{\omega \left(i+k_z^o{\Lambda }^o\right)\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)}{{\Lambda }^o}& -k_x\left(N^{o2}\left(i+k_z^o{\Lambda }^o\right)+\omega \left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\right)& 0& 0\\ 0& {\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega & \frac{i{k}_z^o{\mathrm{e}}^{\frac{z}{{\Lambda }^o}}}{{\rho }_{00}^o}& 0\\ i{k}_x& i{k}_z^o-\frac{1}{{\Lambda }^o}& 0& -i\omega \\ i{k}_x& -\frac{g}{c^{o2}}+i{k}_z^o& -\frac{i\omega }{{\rho }_{00}^o{c}^{o2}}& 0\\ -g{k}_x^2+\frac{\omega \left(i+k_z^a{\Lambda }^a\right)\left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)}{{\Lambda }^a}& -k_x\left(N^{a2}\left(i+k_z^a{\Lambda }^a\right)+\omega \left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\right)& 0& 0\\ 0& {\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega & \frac{i{k}_z^a{\mathrm{e}}^{\frac{z}{{\Lambda }^a}}}{{\rho }_{00}^a}& 0\\ i{k}_x& i{k}_z^a-\frac{1}{{\Lambda }^a}& 0& -i\omega \\ i{k}_x& -\frac{g}{c^{a2}}+i{k}_z^a& -\frac{i\omega }{{\rho }_{00}^a{c}^{a2}}& 0\end{array}\right)=\left(\begin{array}{l}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)$$

(55)

The resulting system is divided into two independent systems of equations that describe the relationships for the upper and lower media. The compatibility condition for each of the systems leads to dispersion relations for the lower one:

$$\begin{array}{c}\frac{\omega }{c^{o2}{\Lambda }^{o2}}\left[\omega \left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\left(-g{k}_x^2{\Lambda }^o+\omega \left(i+k_z^0{\Lambda }^o\right)\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\right)+\right.\\ +{\mathrm{e}}^{\frac{z}{{\Lambda }^o}}{k}_z^o\left(\left(g+i{c}^{o2}{k}_z^{o2}\right)\left(-g{k}_x^2{\Lambda }^o+\omega \left(i+k_z^0{\Lambda }^o\right)\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\right)\right.+\\ +\left.\left.c^{o2}{k}_x^2{\Lambda }^o\left(N^{o2}\left(i{k}_z^o{\Lambda }^o-1\right)+\omega \left(i\nu \left(k_x^2+k_z^{o2}\right)+\omega \right)\right)\right)\right]=0\end{array}$$

(56)

and top liquid:

$$\begin{array}{c}\frac{\omega }{c^{a2}{\Lambda }^{a2}}\left[\omega \left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\left(-g{k}_x^2{\Lambda }^a+\omega \left(i+k_z^a{\Lambda }^a\right)\left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\right)+\right.\\ +{\mathrm{e}}^{\frac{z}{{\Lambda }^a}}{k}_z^a\left(\left(g+i{c}^{a2}{k}_z^{a2}\right)\left(-g{k}_x^2{\Lambda }^a+\omega \left(i+k_z^a{\Lambda }^a\right)\left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\right)\right.+\\ +\left.\left.c^{a2}{k}_x^2{\Lambda }^a\left(N^{a2}\left(i{k}_z^a{\Lambda }^a-1\right)+\omega \left(i\nu \left(k_x^2+k_z^{a2}\right)+\omega \right)\right)\right)\right]=0\end{array}$$

(57)

Let us consider expressions (56) – (57) in dimensionless form. As non-dimensional scales we will choose, as in the previous paragraphs, the natural parameters of the medium: as the time scale, we will choose the inverse buoyancy frequency ${\tau }_b^{o,a}=N^{o,a-1}$, and as the spatial scale, we will choose the viscous wave scale ${\delta }_N^{g\nu o,a}={\left(g{\nu }^{o,a}\right)}^{1/3}{N}^{o,a-1}$:

$$\begin{array}{c}\frac{{\omega }_{*}}{{\epsilon }^{o2}}\left[{\epsilon }^o{\eta }^o{\omega }_{*}\left({\epsilon }^o{\omega }_{*}\left(k_{z*}^o+i{\epsilon }^o\right)\left({\epsilon }^o{k}_{z*}^{o2}-i{\omega }_{*}\right)+{\epsilon }^o{k}_{x*}^4\left(-1+{\epsilon }^{o2}{\omega }_{*}{k}_{z*}^o+i{\epsilon }^{o3}{\omega }_{*}\right)\right.\right.+\\ \left.+k_{x*}^2\left(2{\epsilon }^{o3}{\omega }_{*}{k}_{z*}^{o3}-2i{\epsilon }^{o2}{\omega }_{*}^2{k}_{z*}^o+{\omega }_{*}\left(i+2{\epsilon }^{o3}{\omega }_{*}\right)+k_{z*}^{o2}\left(-{\epsilon }^o+2i{\epsilon }^{o4}{\omega }_{*}\right)\right)\right)+\\ +{\mathrm{e}}^{\frac{z}{{\Lambda }^o}}{k}_{z*}^o\left(i{\epsilon }^{o4}{\omega }_{*}{k}_{x*}^4+\right.{\epsilon }^o{\omega }_{*}\left(k_{z*}^o+i{\epsilon }^o\right)\left({\epsilon }^{o2}{k}_{z*}^o-i{\eta }^o\right)\left(i{\epsilon }^o{k}_{z*}^{o2}+{\omega }_{*}\right)-\\ -\left.\left.k_{x*}^2\left({\eta }^o-2i{\epsilon }^{o4}{\omega }_{*}{k}_{z*}^{o2}+{\epsilon }^{o5}{\omega }_{*}{k}_{z*}^o-{\epsilon }^{o2}{\eta }^o{\omega }_{*}{k}_{z*}^o+{\epsilon }^{o3}\left(1-i{\eta }^o{\omega }_{*}-{\omega }_{*}^2\right)\right)\right)\right]=0\end{array}$$

(58)

$$\begin{array}{c}\frac{{\omega }_{*}}{{\epsilon }^{a2}}\left[{\epsilon }^a{\eta }^a{\omega }_{*}\left({\epsilon }^a{\omega }_{*}\left(k_{z*}^a+i{\epsilon }^a\right)\left({\epsilon }^a{k}_{z*}^{a2}-i{\omega }_{*}\right)+{\epsilon }^a{k}_{x*}^4\left(-1+{\epsilon }^{a2}{\omega }_{*}{k}_{z*}^a+i{\epsilon }^{a3}{\omega }_{*}\right)\right.\right.+\\ \left.+k_{x*}^2\left(2{\epsilon }^{a3}{\omega }_{*}{k}_{z*}^{a3}-2i{\epsilon }^{a2}{\omega }_{*}^2{k}_{z*}^a+{\omega }_{*}\left(i+2{\epsilon }^{a3}{\omega }_{*}\right)+k_{z*}^{a2}\left(-{\epsilon }^a+2i{\epsilon }^{a4}{\omega }_{*}\right)\right)\right)+\\ +{\mathrm{e}}^{\frac{z}{{\Lambda }^a}}{k}_{z*}^a\left(i{\epsilon }^{a4}{\omega }_{*}{k}_{x*}^4+\right.{\epsilon }^a{\omega }_{*}\left(k_{z*}^a+i{\epsilon }^a\right)\left({\epsilon }^{a2}{k}_{z*}^a-i{\eta }^a\right)\left(i{\epsilon }^a{k}_{z*}^{a2}+{\omega }_{*}\right)-\\ -\left.\left.k_{x*}^2\left({\eta }^a-2i{\epsilon }^{a4}{\omega }_{*}{k}_{z*}^{a2}+{\epsilon }^{a5}{\omega }_{*}{k}_{z*}^a-{\epsilon }^{a2}{\eta }^a{\omega }_{*}{k}_{z*}^a+{\epsilon }^{a3}\left(1-i{\eta }^a{\omega }_{*}-{\omega }_{*}^2\right)\right)\right)\right]=0\end{array}$$

(59)

$${\epsilon }^a=N^a\sqrt[3]{\frac{{\nu }^a}{g^2}},{\epsilon }^o=N^o\sqrt[3]{\frac{{\nu }^o}{g^2}},{\eta }^a=\frac{N^a{\nu }^a}{c^{a2}},{\eta }^o=\frac{N^o{\nu }^o}{c^{o2}}.$$

Expressions (58) – (59) are reduced to the dispersion relations in an incompressible fluid when passing to the limit $c^{o,a}\to \infty$ (${\eta }^{o,a}\to 0$):

$${\omega }_{*}\left(i{\epsilon }^o\left(k_{x*}^2+k_{z*}^{o2}\right)+{\omega }_{*}\right)\left({\epsilon }^o{\omega }_{*}\left(k_{z*}^o+i{\epsilon }^o\right)\left({\epsilon }^o{k}_{z*}^{o2}-i{\omega }_{*}\right)+k_{x*}^2\left(-1+{\epsilon }^{o2}{\omega }_{*}{k}_{z*}^o+i{\epsilon }^{o3}{\omega }_{*}\right)\right)=0$$

(60)

$${\omega }_{*}\left(i{\epsilon }^a\left(k_{x*}^2+k_{z*}^{a2}\right)+{\omega }_{*}\right)\left({\epsilon }^a{\omega }_{*}\left(k_{z*}^a+i{\epsilon }^a\right)\left({\epsilon }^a{k}_{z*}^{a2}-i{\omega }_{*}\right)+k_{x*}^2\left(-1+{\epsilon }^{a2}{\omega }_{*}{k}_{z*}^a+i{\epsilon }^{a3}{\omega }_{*}\right)\right)=0$$

(61)

The small parameter ${\eta }^{o,a}$for liquids with the parameters of water and air turns out to be significantly smaller than the small parameter ${\epsilon }^{o,a}$. Approximate solutions of dispersion relations (58) – (59) relatively $k_{*z}^{o,a}$have the form:

$$k_{*z}^{o,a}=k_{0*z}^{o,a}+\eta k_{1*z}^{o,a}$$

(62)

In solution $k_{0*z}^{o,a}$ (62) takes one of the following values:

$$k_{0*z}^{o,a}=0;$$

(63)

$$k_{0*z}^{o,a}=-\frac{i{\epsilon }^{o,a}}{4}-\frac{1}{2}\sqrt{-\frac{{\epsilon }^{o,a}}{4}-\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{{\epsilon }^{o,a}}+\theta }\pm \frac{1}{2}\sqrt{-\frac{{\epsilon }^{o,a2}}{2}-\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{{\epsilon }^{o,a}}-\theta -\frac{i{\epsilon }^{o,a3}-8i\left({\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}\right)+4i\left(2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}\right)}{4\sqrt{-\frac{{\epsilon }^{o,a2}}{4}+\theta -\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{{\epsilon }^{o,a}}}}}$$

(64)

$$k_{0*z}^{o,a}=-\frac{i{\epsilon }^{o,a}}{4}+\frac{1}{2}\sqrt{-\frac{{\epsilon }^{o,a}}{4}-\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{{\epsilon }^{o,a}}+\theta }\pm \frac{1}{2}\sqrt{-\frac{{\epsilon }^{o,a2}}{2}-\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{{\epsilon }^{o,a}}-\theta +\frac{i{\epsilon }^{o,a3}-8i\left({\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}\right)+4i\left(2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}\right)}{4\sqrt{-\frac{{\epsilon }^{o,a2}}{4}-\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{{\epsilon }^{o,a}}+\theta }}};$$

(65)

$$\theta =\frac{\left(i+\sqrt{3}\right){\left(\alpha +\sqrt{{\alpha }^2-4{\beta }^3}\right)}^{1/3}}{6\cdot 2^{1/3}{\epsilon }^{o,a}{\omega }_{*}}+\frac{2{\epsilon }^{o,a}{k}_{x*}^2-i{\omega }_{*}}{3{\epsilon }^{o,a}}+\frac{\left(i-\sqrt{3}\right)2^{1/3}\beta }{3{\epsilon }^{o,a}{\left(\alpha +\sqrt{{\alpha }^2-4{\beta }^3}\right)}^{1/3}}$$

$$\beta =\left(-16{\epsilon }^{o,a2}{k}_{x*}^4{\omega }_{*}+{\omega }_{*}^2\left(3i{\epsilon }^{o,a}+{\omega }_{*}\right)+k_{x*}^2\left(-3{\epsilon }^{o,a4}{\omega }_{*}+4i{\epsilon }^{o,a}\left(-3+4{\omega }_{*}^2\right)\right)\right)$$

$$\alpha ={\omega }_{*}^2\left(128i{\epsilon }^{o,a3}{\omega }_{*}{k}_{x*}^6+2{\omega }_{*}^3\left(-9i{\epsilon }^{o,a3}+{\omega }_{*}\right)+12{\epsilon }^{o,a2}{k}_{x*}^4\left(-12+3i{\epsilon }^{o,a3}{\omega }_{*}+16{\omega }_{*}^2\right)+3k_{x*}^2\left(9{\epsilon }^{o,a4}\left(-1+2{\omega }_{*}^2\right)-4i{\epsilon }^{o,a}{\omega }_{*}\left(-6+5{\omega }_{*}^2\right)\right)\right)$$

and$k_{1*z}^{o,a}$ takes the corresponding (63) – (65) values:

$$k_{1*z}^{o,a}=\frac{\left(k_{0z*}^{o,a}+{\epsilon }^{o,a}{\mathrm{e}}^{-\frac{z}{{\Lambda }^{o,a}}}\left({\epsilon }^{o,a}{\omega }_{*}\left(k_{x*}^2+k_{0z*}^{o,a2}-i{\omega }_{*}\right)\right)\right)\left({\epsilon }^{o,a}{\omega }_{*}\left(k_{0z*}^{o,a}+i{\epsilon }^{o,a}\right)\left({\epsilon }^{o,a}{k}_{0z*}^{o,a2}-i{\omega }_{*}\right)+k_{x*}^2\left(-1+{\epsilon }^{o,a2}{k}_{0z*}^{o,a}{\omega }_{*}+i{\epsilon }^{o,a3}{\omega }_{*}\right)\right)}{{\epsilon }^{o,a3}\left(-i{\epsilon }^{o,a}{\omega }_{*}{k}_{x*}^4+{\omega }_{*}{k}_{0z*}^{o,a}\left(-5i{\epsilon }^{o,a}{k}_{0z*}^{o,a3}+4{\epsilon }^{o,a2}{k}_{0z*}^2-3{\omega }_{*}{k}_{0z*}^{o,a}-2i{\epsilon }^{o,a}{\omega }_{*}\right)+k_{x*}^2\left(1-6i{\epsilon }^{o,a}{\omega }_{*}{k}_{0z*}^{o,a2}+2{\epsilon }^{o,a2}{\omega }_{*}{k}_{0z*}^{o,a}-{\omega }_{*}^2\right)\right)}$$

(66)

Additional conditions for physical implementation are imposed on solutions (63) – (66):

$$\mathrm{Im}\left(k_{*x}\right)>0,\mathrm{Im}\left(k_{*z}^o\right)<0,\mathrm{Im}\left(k_{*z}^a\right)>0$$

(67)

Taking (67) into account, solution (63) turns out to be physically unrealizable in both media. Solution (64) describes a regular solution with respect to a small parameter ${\epsilon }^{o,a}$ and the corresponding wave component of a periodic flow. Solution (65) describes a singular solution with respect to a small parameter ${\epsilon }^{o,a}$and correspond to the ligament component of the periodic flow. To distinguish the roots, we introduce a redesignation for singular solutions $k_{*l}^{o,a}$. Mathematically, solutions corresponding to the wave component are determined by the condition:

$$\left|\mathrm{Re}\left(k_{*z}^{o,a}\right)\right|\gg \left|\mathrm{Im}\left(k_{*z}^{o,a}\right)\right|$$

(68)

and solutions corresponding to the ligament component are determined by the mathematical condition:

$$\left|\mathrm{Re}\left(k_{*l}^{o,a}\right)\right|\sim \left|\mathrm{Im}\left(k_{*l}^{o,a}\right)\right|$$

(69)

Taking into account the ligament components, the form of the complete solution (54) will be rewritten as:

$$\left(\begin{array}{l}{u}^{o,a}\\ w^{o,a}\\ {\tilde{P}}^{o,a}\\ {\tilde{\rho }}^{o,a}\\ \zeta \end{array}\right)=\left(\begin{array}{c}{U}_m^{o,a}\left(\exp \left(i{k}_z^{o,a}z\right)+\Theta \exp \left(i{k}_l^{o,a}z\right)\right)\\ W_m^{o,a}\left(\exp \left(i{k}_z^{o,a}z\right)+\Theta \exp \left(i{k}_l^{o,a}z\right)\right)\\ P_m^{o,a}\left(\exp \left(i{k}_z^{o,a}z\right)+\Theta \exp \left(i{k}_l^{o,a}z\right)\right)\\ {{\rm P}}_m^{o,a}\left(\exp \left(i{k}_z^{o,a}z\right)+\Theta \exp \left(i{k}_l^{o,a}z\right)\right)\\ A_m\end{array}\right)\exp \left(i{k}_{x}x-i\omega t\right)$$

(70)

Substituting the form of solution (70) into the boundary conditions (53), we obtain dispersion relations connecting the components of the wave vector $k_x$ with the frequency of wave motion $\omega$. Substituting approximate solutions (64), (66) and (65), (66) into the resulting relation we obtain a dispersion equation. Restrictions (67) are imposed on the solution, thus physically realizable roots are selected. The resulting expressions are cumbersome and difficult to analyze. Let's consider some limiting cases.

Let us consider the behavior of oscillations far from the interface between the media. In this case, we will assume that $\left|z\right|\gg 1$. For the lower liquid, for the upper -. Thus, for the ocean in the dispersion relation (56), we can neglect the second term and remain:

$$\frac{{\omega }^2\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\left(-g{k}_x^2{\Lambda }^o+\omega \left(i+k_z^0{\Lambda }^o\right)\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\right)}{c^{o2}{\Lambda }^{o2}}=0$$

(71)

Or in disdimention form:

$${\omega }_{*}^2{\eta }^o\left({\omega }_{*}\left(k_{z*}^o+i{\epsilon }^o\right)\left({\epsilon }^o{k}_{z*}^{o2}-i{\omega }_{*}\right)+k_{x*}^4\left(-1+{\epsilon }^{o2}{\omega }_{*}{k}_{z*}^o+i{\epsilon }^{o3}{\omega }_{*}\right)\right.=0$$

(72)

The solutions to expression (71) (or (72)) are found exactly, but due to their cumbersomeness they are not given here.

For the atmosphere in the dispersion relation (57), from similar reasoning we can neglect the first term and obtain the dispersion relation far from the interface:

$$\frac{\omega k_z^a\left(\left(g+i{c}^{a2}{k}_z^{a2}\right)\left(-g{k}_x^2{\Lambda }^a+\omega \left(i+k_z^a{\Lambda }^a\right)\left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\right)+c^{a2}{k}_x^2{\Lambda }^a\left(N^{a2}\left(i{k}_z^a{\Lambda }^a-1\right)+\omega \left(i\nu \left(k_x^2+k_z^{a2}\right)+\omega \right)\right)\right)}{c^{a2}{\Lambda }^{a2}}=0$$

(73)

Or in disdimention form:

$$\begin{array}{c}\frac{{\omega }_{*}^2}{{\epsilon }^a}{\eta }^a\left({\epsilon }^a{\omega }_{*}\left(k_{z*}^a+i{\epsilon }^a\right)\left({\epsilon }^a{k}_{z*}^{a2}-i{\omega }_{*}\right)+{\epsilon }^a{k}_{x*}^4\left(-1+{\epsilon }^{a2}{\omega }_{*}{k}_{z*}^a+i{\epsilon }^{a3}{\omega }_{*}\right)+\right.\\ \left.+k_{x*}^2\left(2{\epsilon }^{a3}{\omega }_{*}{k}_{z*}^{a3}-2i{\epsilon }^{a2}{\omega }_{*}^2{k}_{z*}^a+{\omega }_{*}\left(i+2{\epsilon }^{a3}{\omega }_{*}\right)+k_{z*}^{a2}\left(-{\epsilon }^a+2i{\epsilon }^{a4}{\omega }_{*}\right)\right)\right)=0\end{array}$$

(74)

The solutions to expression (73) (or (74)) also due to their cumbersomeness are not given here. For waves near the surface, we can assume that dispersion relations (56) – (57) will be simplified:

$$\begin{array}{c}\frac{\omega }{c^{o2}{\Lambda }^{o2}}\left[\omega \left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\left(-g{k}_x^2{\Lambda }^o+\omega \left(i+k_z^0{\Lambda }^o\right)\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\right)+\right.\\ +k_z^o\left(\left(g+i{c}^{o2}{k}_z^{o2}\right)\left(-g{k}_x^2{\Lambda }^o+\omega \left(i+k_z^0{\Lambda }^o\right)\left({\nu }^o\left(k_x^2+k_z^{o2}\right)-i\omega \right)\right)\right.+\\ +\left.\left.c^{o2}{k}_x^2{\Lambda }^o\left(N^{o2}\left(i{k}_z^o{\Lambda }^o-1\right)+\omega \left(i\nu \left(k_x^2+k_z^{o2}\right)+\omega \right)\right)\right)\right]=0\end{array}$$

(75)

$$\begin{array}{c}\frac{\omega }{c^{a2}{\Lambda }^{a2}}\left[\omega \left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\left(-g{k}_x^2{\Lambda }^a+\omega \left(i+k_z^a{\Lambda }^a\right)\left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\right)+\right.\\ +k_z^a\left(\left(g+i{c}^{a2}{k}_z^{a2}\right)\left(-g{k}_x^2{\Lambda }^a+\omega \left(i+k_z^a{\Lambda }^a\right)\left({\nu }^a\left(k_x^2+k_z^{a2}\right)-i\omega \right)\right)\right.+\\ +\left.\left.c^{a2}{k}_x^2{\Lambda }^a\left(N^{a2}\left(i{k}_z^a{\Lambda }^a-1\right)+\omega \left(i\nu \left(k_x^2+k_z^{a2}\right)+\omega \right)\right)\right)\right]=0\end{array}$$

(76)

Or in disdimention form:

$$\begin{array}{c}\frac{{\omega }_{*}}{{\epsilon }^{o2}}\left[{\epsilon }^o{\eta }^o{\omega }_{*}\left({\epsilon }^o{\omega }_{*}\left(k_{z*}^o+i{\epsilon }^o\right)\left({\epsilon }^o{k}_{z*}^{o2}-i{\omega }_{*}\right)+{\epsilon }^o{k}_{x*}^4\left(-1+{\epsilon }^{o2}{\omega }_{*}{k}_{z*}^o+i{\epsilon }^{o3}{\omega }_{*}\right)\right.\right.+\\ \left.+k_{x*}^2\left(2{\epsilon }^{o3}{\omega }_{*}{k}_{z*}^{o3}-2i{\epsilon }^{o2}{\omega }_{*}^2{k}_{z*}^o+{\omega }_{*}\left(i+2{\epsilon }^{o3}{\omega }_{*}\right)+k_{z*}^{o2}\left(-{\epsilon }^o+2i{\epsilon }^{o4}{\omega }_{*}\right)\right)\right)+\\ +k_{z*}^o\left(i{\epsilon }^{o4}{\omega }_{*}{k}_{x*}^4+\right.{\epsilon }^o{\omega }_{*}\left(k_{z*}^o+i{\epsilon }^o\right)\left({\epsilon }^{o2}{k}_{z*}^o-i{\eta }^o\right)\left(i{\epsilon }^o{k}_{z*}^{o2}+{\omega }_{*}\right)-\\ -\left.\left.k_{x*}^2\left({\eta }^o-2i{\epsilon }^{o4}{\omega }_{*}{k}_{z*}^{o2}+{\epsilon }^{o5}{\omega }_{*}{k}_{z*}^o-{\epsilon }^{o2}{\eta }^o{\omega }_{*}{k}_{z*}^o+{\epsilon }^{o3}\left(1-i{\eta }^o{\omega }_{*}-{\omega }_{*}^2\right)\right)\right)\right]=0\end{array}$$

(77)

$$\begin{array}{c}\frac{{\omega }_{*}}{{\epsilon }^{a2}}\left[{\epsilon }^a{\eta }^a{\omega }_{*}\left({\epsilon }^a{\omega }_{*}\left(k_{z*}^a+i{\epsilon }^a\right)\left({\epsilon }^a{k}_{z*}^{a2}-i{\omega }_{*}\right)+{\epsilon }^a{k}_{x*}^4\left(-1+{\epsilon }^{a2}{\omega }_{*}{k}_{z*}^a+i{\epsilon }^{a3}{\omega }_{*}\right)\right.\right.+\\ \left.+k_{x*}^2\left(2{\epsilon }^{a3}{\omega }_{*}{k}_{z*}^{a3}-2i{\epsilon }^{a2}{\omega }_{*}^2{k}_{z*}^a+{\omega }_{*}\left(i+2{\epsilon }^{a3}{\omega }_{*}\right)+k_{z*}^{a2}\left(-{\epsilon }^a+2i{\epsilon }^{a4}{\omega }_{*}\right)\right)\right)+\\ +k_{z*}^a\left(i{\epsilon }^{a4}{\omega }_{*}{k}_{x*}^4+\right.{\epsilon }^a{\omega }_{*}\left(k_{z*}^a+i{\epsilon }^a\right)\left({\epsilon }^{a2}{k}_{z*}^a-i{\eta }^a\right)\left(i{\epsilon }^a{k}_{z*}^{a2}+{\omega }_{*}\right)-\\ -\left.\left.k_{x*}^2\left({\eta }^a-2i{\epsilon }^{a4}{\omega }_{*}{k}_{z*}^{a2}+{\epsilon }^{a5}{\omega }_{*}{k}_{z*}^a-{\epsilon }^{a2}{\eta }^a{\omega }_{*}{k}_{z*}^a+{\epsilon }^{a3}\left(1-i{\eta }^a{\omega }_{*}-{\omega }_{*}^2\right)\right)\right)\right]=0\end{array}$$

(78)

But despite their simpler appearance, the roots of expressions (77) – (78), as well as complete expressions, can only be found asymptotically or numerically.

7. Discussion

The expressive properties of periodic flows in fluids - the regularity of wave displacements of the free surface of the liquid, the high speed of propagation of sound vibrations, the clarity of the pattern of beams of periodic internal waves - formed the basis for the generally accepted classification of waves and predetermined the rules for constructing mathematical models of the phenomenon. To describe each wave process in a linear [1, 2, 16, 21] or nonlinear approximation [41], its own system of equations was developed based on the system of fundamental equations of mechanics of fluids and gases [1, 2, 4, 16] , and general physical considerations [31,50].

Under natural conditions, sharp disturbances lead to the formation of several types of waves, which propagate with their own phase and group velocities and differ in attenuation laws. The parameters of wave processes - periods, wavelength, group and phase propagation velocities are described by real numbers. The mathematical description of periodic flows is carried out in the algebra of complex numbers. The use of wave representations by exponential functions of complex frequency and complex wave vector allows us to construct dispersion relations [1, 2] and evaluate the stability of the flows under study [32,33].

Taking into account the special physical properties of the wave frequency - a measure of the energy of periodic motion, in this work, as in [38, 40, 42, 44], the wave frequency is assumed to be real, and the wave number is taken to be complex. In this approximation, the degree of the dispersion relation corresponds to the order of the system of differential equations. Solutions of the system of governing equations, constructed by methods of singular perturbation theory taking into account the type of small parameter of the process under study, contain two types of solutions. The real part of some wave numbers is large, and the imaginary part is small. Others have real and imaginary parts of the same order. Accordingly, some of the solutions, including solutions with small values of the imaginary part of the wave vectors, contain functions that are regular in the small parameter and describe waves. For each type of wave, its own dispersion equation is constructed.

Another part of the solutions with large values of the imaginary parts of the wave vector determines the ligaments, which correspond to thin high-gradient fibers and interfaces in the thickness of the stratified liquid [30, 38]. From the analysis it follows that specific ligaments accompany all types of waves - surface, internal, and acoustic. Taking into account the influence of ligaments made it possible to precalculate the parameters of reflected and leaking waves when reflecting beams of internal waves from the critical level separating a medium with a high frequency of buoyancy from a layer with a low frequency not exceeding the frequency of the wave [5], consistent with the data of later experiments [51].

In theory, the number of ligaments accompanying the wave is determined by the degree of completeness of taking into account the factors influencing the density and the dimension of the problem space. The minimum number - two ligaments - accompany two-dimensional waves in a medium with one dissipative parameter (kinematic viscosity). Their thickness is determined by the scale of the periodic Stokes flow ${\delta }_{\omega }^{\nu }=\sqrt{\nu /\omega }$ [13]. Taking into account the three-dimensionality of space, the effects of thermal diffusivity and diffusion leads to an increase in the number of ligaments with different properties [40]. The effects of nonlinear interaction between ligaments can increase the mutual influence of waves of different types [42].

The developed methodology for constructing complete solutions makes it possible to describe not only the wave component of a periodic flow, but also the fine structure, manifested in the form of ligaments - thin jets accompanying the wave motion. The parameters of the observed phenomena in the process of propagation of periodic disturbances in liquids and gases, which are determined by the properties of the medium, determine the requirements for the experimental methodology and the resolution (spatial and temporal) of the equipment for observing the complete picture of flows.